Performance/bin/LU.c

0001 #include <math.h>
0002 #include "LU.h"
0003
0004 double LU_num_flops(int N) {
0005   /* rougly 2/3*N^3 */
0006
0007   double Nd = (double)N;
0008
0009   return (2.0 * Nd * Nd * Nd / 3.0);
0010 }
0011
0012 void LU_copy_matrix(int M, int N, double **lu, double **A) {
0013   int i;
0014   int j;
0015
0016   for (i = 0; i < M; i++)
0017     for (j = 0; j < N; j++)
0018       lu[i][j] = A[i][j];
0019 }
0020
0021 int LU_factor(int M, int N, double **A, int *pivot) {
0022   int minMN = M < N ? M : N;
0023   int j = 0;
0024
0025   for (j = 0; j < minMN; j++) {
0026     /* find pivot in column j and  test for singularity. */
0027
0028     int jp = j;
0029     int i;
0030
0031     double t = fabs(A[j][j]);
0032     for (i = j + 1; i < M; i++) {
0033       double ab = fabs(A[i][j]);
0034       if (ab > t) {
0035         jp = i;
0036         t = ab;
0037       }
0038     }
0039
0040     pivot[j] = jp;
0041
0042     /* jp now has the index of maximum element  */
0043     /* of column j, below the diagonal          */
0044
0045     if (A[jp][j] == 0)
0046       return 1; /* factorization failed because of zero pivot */
0047
0048     if (jp != j) {
0049       /* swap rows j and jp */
0050       double *tA = A[j];
0051       A[j] = A[jp];
0052       A[jp] = tA;
0053     }
0054
0055     if (j < M - 1) /* compute elements j+1:M of jth column  */
0056     {
0057       /* note A(j,j), was A(jp,p) previously which was */
0058       /* guarranteed not to be zero (Label #1)         */
0059
0060       double recp = 1.0 / A[j][j];
0061       int k;
0062       for (k = j + 1; k < M; k++)
0063         A[k][j] *= recp;
0064     }
0065
0066     if (j < minMN - 1) {
0067       /* rank-1 update to trailing submatrix:   E = E - x*y; */
0068       /* E is the region A(j+1:M, j+1:N) */
0069       /* x is the column vector A(j+1:M,j) */
0070       /* y is row vector A(j,j+1:N)        */
0071
0072       int ii;
0073       for (ii = j + 1; ii < M; ii++) {
0074         double *Aii = A[ii];
0075         double *Aj = A[j];
0076         double AiiJ = Aii[j];
0077         int jj;
0078         for (jj = j + 1; jj < N; jj++)
0079           Aii[jj] -= AiiJ * Aj[jj];
0080       }
0081     }
0082   }
0083
0084   return 0;
0085 }