1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      72
      73
      74
      75
      76
      77
      78
      79
      80
      81
      82
      83
      84
      85
      86
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
     119
     120
     121
     122
     123
     124
     125
     126
     127
     128
     129
     130
#include <ulmblas/level1.h>
#include <ulmblas/level2.h>

//-- LU factorization ----------------------------------------------------------

void
lu_unblk_var1(int m, int n, double *A, int incRowA, int incColA)
{
    int min_mn = (m<n) ? m : n;
    int max_mn = (m>n) ? m : n;
    int j;

    for (j=1; j<max_mn; ++j) {
        int k = (min_mn < j) ? min_mn : j;

        if (j<n) {
            dtrlsv(k, 1,
                   A, incRowA, incColA,
                   &A[j*incColA], incRowA);
        }
        if (j<m) {
            dtrlsv(k, 0,
                   A, incColA, incRowA,
                   &A[j*incRowA], incColA);
        }
        if (j<m && j<n) {
            A[j*(incRowA+incColA)] -= ddot(j,
                                           &A[j*incRowA], incColA,
                                           &A[j*incColA], incRowA);
        }
    }
}

void
lu_unblk_var2(int m, int n, double *A, int incRowA, int incColA)
{
    int min_mn = (m<n) ? m : n;
    int j;

    for (j=1; j<m; ++j) {
        int k = (min_mn < j) ? min_mn : j;

        if (j<m) {
            dtrlsv(k, 0,
                   A, incColA, incRowA,
                   &A[j*incRowA], incColA);
        }
        if (j<n) {
            dgemv(n-j, j,
                  -1.0,
                  &A[j*incColA], incColA, incRowA,
                  &A[j*incRowA], incColA,
                  1.0,
                  &A[j*(incRowA+incColA)], incColA);
        }
    }
}

void
lu_unblk_var3(int m, int n, double *A, int incRowA, int incColA)
{
    int min_mn = (m<n) ? m : n;
    int j;

    for (j=0; j<n; ++j) {
        int k = (min_mn < j) ? min_mn : j;

        dtrlsv(k, 1,
               A, incRowA, incColA,
               &A[j*incColA], incRowA);
        if (j<m) {
            dgemv(m-j, j,
                  -1.0,
                  &A[j*incRowA], incRowA, incColA,
                  &A[j*incColA], incRowA,
                  1.0,
                  &A[j*(incRowA+incColA)], incRowA);
            if (j+1<m) {
                dscal(m-j-1,
                      1.0/A[j*(incRowA+incColA)],
                      &A[(j+1)*incRowA+j*incColA], incRowA);
            }
        }
    }
}

void
lu_unblk_var4(int m, int n, double *A, int incRowA, int incColA)
{
    int min_mn = (m<n) ? m : n;
    int j;

    for (j=0; j<min_mn; ++j) {
        dgemv(n-j, j,
              -1.0,
              &A[j*incColA], incColA, incRowA,
              &A[j*incRowA], incColA,
              1.0,
              &A[j*(incRowA+incColA)], incColA);
        if (j+1<m) {
            dgemv(m-j-1, j,
                  -1.0,
                  &A[(j+1)*incRowA], incRowA, incColA,
                  &A[j*incColA], incRowA,
                  1.0,
                  &A[(j+1)*incRowA+j*incColA], incRowA);
            dscal(m-j-1,
                  1.0/A[j*(incRowA+incColA)],
                  &A[(j+1)*incRowA+j*incColA], incRowA);
        }
    }
}

void
lu_unblk_var5(int m, int n, double *A, int incRowA, int incColA)
{
    int min_mn = (m<n) ? m : n;
    int j;

    for (j=0; j<min_mn; ++j) {
        dscal(m-j-1,
              1.0/A[j*(incRowA+incColA)],
              &A[(j+1)*incRowA+j*incColA], incRowA);
        dger(m-j-1, n-j-1,
             -1.0,
             &A[(j+1)*incRowA+j*incColA], incRowA,
             &A[j*incRowA+(j+1)*incColA], incColA,
             &A[(j+1)*(incRowA+incColA)], incRowA, incColA);
    }
}