00001 #include "cs.h"
00002
00003 csn *cs_lu (const cs *A, const css *S, double tol)
00004 {
00005 cs *L, *U ;
00006 csn *N ;
00007 double pivot, *Lx, *Ux, *x, a, t ;
00008 int *Lp, *Li, *Up, *Ui, *pinv, *xi, *q, n, ipiv, k, top, p, i, col, lnz,unz;
00009 if (!CS_CSC (A) || !S) return (NULL) ;
00010 n = A->n ;
00011 q = S->q ; lnz = S->lnz ; unz = S->unz ;
00012 x = cs_malloc (n, sizeof (double)) ;
00013 xi = cs_malloc (2*n, sizeof (int)) ;
00014 N = cs_calloc (1, sizeof (csn)) ;
00015 if (!x || !xi || !N) return (cs_ndone (N, NULL, xi, x, 0)) ;
00016 N->L = L = cs_spalloc (n, n, lnz, 1, 0) ;
00017 N->U = U = cs_spalloc (n, n, unz, 1, 0) ;
00018 N->pinv = pinv = cs_malloc (n, sizeof (int)) ;
00019 if (!L || !U || !pinv) return (cs_ndone (N, NULL, xi, x, 0)) ;
00020 Lp = L->p ; Up = U->p ;
00021 for (i = 0 ; i < n ; i++) x [i] = 0 ;
00022 for (i = 0 ; i < n ; i++) pinv [i] = -1 ;
00023 for (k = 0 ; k <= n ; k++) Lp [k] = 0 ;
00024 lnz = unz = 0 ;
00025 for (k = 0 ; k < n ; k++)
00026 {
00027
00028 Lp [k] = lnz ;
00029 Up [k] = unz ;
00030 if ((lnz + n > L->nzmax && !cs_sprealloc (L, 2*L->nzmax + n)) ||
00031 (unz + n > U->nzmax && !cs_sprealloc (U, 2*U->nzmax + n)))
00032 {
00033 return (cs_ndone (N, NULL, xi, x, 0)) ;
00034 }
00035 Li = L->i ; Lx = L->x ; Ui = U->i ; Ux = U->x ;
00036 col = q ? (q [k]) : k ;
00037 top = cs_spsolve (L, A, col, xi, x, pinv, 1) ;
00038
00039 ipiv = -1 ;
00040 a = -1 ;
00041 for (p = top ; p < n ; p++)
00042 {
00043 i = xi [p] ;
00044 if (pinv [i] < 0)
00045 {
00046 if ((t = fabs (x [i])) > a)
00047 {
00048 a = t ;
00049 ipiv = i ;
00050 }
00051 }
00052 else
00053 {
00054 Ui [unz] = pinv [i] ;
00055 Ux [unz++] = x [i] ;
00056 }
00057 }
00058 if (ipiv == -1 || a <= 0) return (cs_ndone (N, NULL, xi, x, 0)) ;
00059 if (pinv [col] < 0 && fabs (x [col]) >= a*tol) ipiv = col ;
00060
00061 pivot = x [ipiv] ;
00062 Ui [unz] = k ;
00063 Ux [unz++] = pivot ;
00064 pinv [ipiv] = k ;
00065 Li [lnz] = ipiv ;
00066 Lx [lnz++] = 1 ;
00067 for (p = top ; p < n ; p++)
00068 {
00069 i = xi [p] ;
00070 if (pinv [i] < 0)
00071 {
00072 Li [lnz] = i ;
00073 Lx [lnz++] = x [i] / pivot ;
00074 }
00075 x [i] = 0 ;
00076 }
00077 }
00078
00079 Lp [n] = lnz ;
00080 Up [n] = unz ;
00081 Li = L->i ;
00082 for (p = 0 ; p < lnz ; p++) Li [p] = pinv [Li [p]] ;
00083 cs_sprealloc (L, 0) ;
00084 cs_sprealloc (U, 0) ;
00085 return (cs_ndone (N, NULL, xi, x, 1)) ;
00086 }