Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
lu.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 
8 #include "fasp.h"
9 #include "fasp_functs.h"
10 
11 /*---------------------------------*/
12 /*-- Public Functions --*/
13 /*---------------------------------*/
14 
47  INT pivot[],
48  const INT n)
49 {
50  INT i, j, k;
51  REAL *p_k=NULL, *p_row=NULL, *p_col=NULL;
52  REAL max;
53 
54  /* For each row and column, k = 0, ..., n-1, */
55  for (k = 0, p_k = A; k < n; p_k += n, k++) {
56 
57  // find the pivot row
58  pivot[k] = k;
59  max = fabs( *(p_k + k) );
60  for (j = k + 1, p_row = p_k + n; j < n; ++j, p_row += n) {
61  if ( max < fabs(*(p_row + k)) ) {
62  max = fabs(*(p_row + k));
63  pivot[k] = j;
64  p_col = p_row;
65  }
66  }
67 
68  // if the pivot row differs from the current row, interchange the two rows.
69  if (pivot[k] != k)
70  for (j = 0; j < n; ++j) {
71  max = *(p_k + j);
72  *(p_k + j) = *(p_col + j);
73  *(p_col + j) = max;
74  }
75 
76  // if the matrix is singular, return error
77  if ( fabs( *(p_k + k) ) < SMALLREAL ) return -1;
78 
79  // otherwise find the lower triangular matrix elements for column k.
80  for (i = k+1, p_row = p_k + n; i < n; p_row += n, ++i) {
81  *(p_row + k) /= *(p_k + k);
82  }
83 
84  // update remaining matrix
85  for (i = k+1, p_row = p_k + n; i < n; p_row += n, ++i)
86  for (j = k+1; j < n; ++j)
87  *(p_row + j) -= *(p_row + k) * *(p_k + j);
88 
89  }
90 
91  return FASP_SUCCESS;
92 }
93 
118  REAL b[],
119  INT pivot[],
120  REAL x[],
121  const INT n)
122 {
123  INT i, k;
124  REAL *p_k;
125  REAL dum;
126 
127  /* solve Ly = b */
128  for (k = 0, p_k = A; k < n; p_k += n, k++) {
129  if (pivot[k] != k) {dum = b[k]; b[k] = b[pivot[k]]; b[pivot[k]] = dum; }
130  x[k] = b[k];
131  for (i = 0; i < k; ++i) x[k] -= x[i] * *(p_k + i);
132  }
133 
134  /* solve Ux = y */
135  for (k = n-1, p_k = A + n*(n-1); k >= 0; k--, p_k -= n) {
136  if (pivot[k] != k) {dum = b[k]; b[k] = b[pivot[k]]; b[pivot[k]] = dum; }
137  for (i = k + 1; i < n; ++i) x[k] -= x[i] * *(p_k + i);
138  if (*(p_k + k) == 0.0) return -1;
139  x[k] /= *(p_k + k);
140  }
141 
142  return FASP_SUCCESS;
143 }
144 
145 /*---------------------------------*/
146 /*-- End of File --*/
147 /*---------------------------------*/
148 
149 /*
150 
151  //A simple test example can be written as the following
152  INT main (INT argc, const char * argv[])
153  {
154  REAL A[3][3] = {{0.0, 1.0, 4.0},
155  {4.0, 1.0, 0.0},
156  {1.0, 4.0, 1.0}};
157 
158  REAL b[3] = {1, 1, 1}, x[3];
159 
160  INT pivot[3];
161 
162  INT ret, i, j;
163 
164  ret = lu_decomp(&A[0][0], pivot, 3); // LU decomposition
165 
166  ret = lu_solve(&A[0][0], b, pivot, x, 3); // Solve decomposed Ax=b
167 
168  return 1;
169  }
170 
171 */
#define REAL
Definition: fasp.h:67
#define INT
Definition: fasp.h:64
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
Main header file for FASP.
SHORT fasp_smat_lu_decomp(REAL *A, INT pivot[], const INT n)
LU decomposition of A usind Doolittle's method.
Definition: lu.c:46
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
SHORT fasp_smat_lu_solve(REAL *A, REAL b[], INT pivot[], REAL x[], const INT n)
Solving Ax=b using LU decomposition.
Definition: lu.c:117
#define SMALLREAL
Definition: fasp_const.h:238