00001 #include "cs.h"
00002
00003 csn *cs_qr (const cs *A, const css *S)
00004 {
00005 double *Rx, *Vx, *Ax, *x, *Beta ;
00006 int i, k, p, m, n, vnz, p1, top, m2, len, col, rnz, *s, *leftmost, *Ap, *Ai,
00007 *parent, *Rp, *Ri, *Vp, *Vi, *w, *pinv, *q ;
00008 cs *R, *V ;
00009 csn *N ;
00010 if (!CS_CSC (A) || !S) return (NULL) ;
00011 m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
00012 q = S->q ; parent = S->parent ; pinv = S->pinv ; m2 = S->m2 ;
00013 vnz = S->lnz ; rnz = S->unz ; leftmost = S->leftmost ;
00014 w = cs_malloc (m2+n, sizeof (int)) ;
00015 x = cs_malloc (m2, sizeof (double)) ;
00016 N = cs_calloc (1, sizeof (csn)) ;
00017 if (!w || !x || !N) return (cs_ndone (N, NULL, w, x, 0)) ;
00018 s = w + m2 ;
00019 for (k = 0 ; k < m2 ; k++) x [k] = 0 ;
00020 N->L = V = cs_spalloc (m2, n, vnz, 1, 0) ;
00021 N->U = R = cs_spalloc (m2, n, rnz, 1, 0) ;
00022 N->B = Beta = cs_malloc (n, sizeof (double)) ;
00023 if (!R || !V || !Beta) return (cs_ndone (N, NULL, w, x, 0)) ;
00024 Rp = R->p ; Ri = R->i ; Rx = R->x ;
00025 Vp = V->p ; Vi = V->i ; Vx = V->x ;
00026 for (i = 0 ; i < m2 ; i++) w [i] = -1 ;
00027 rnz = 0 ; vnz = 0 ;
00028 for (k = 0 ; k < n ; k++)
00029 {
00030 Rp [k] = rnz ;
00031 Vp [k] = p1 = vnz ;
00032 w [k] = k ;
00033 Vi [vnz++] = k ;
00034 top = n ;
00035 col = q ? q [k] : k ;
00036 for (p = Ap [col] ; p < Ap [col+1] ; p++)
00037 {
00038 i = leftmost [Ai [p]] ;
00039 for (len = 0 ; w [i] != k ; i = parent [i])
00040 {
00041 s [len++] = i ;
00042 w [i] = k ;
00043 }
00044 while (len > 0) s [--top] = s [--len] ;
00045 i = pinv [Ai [p]] ;
00046 x [i] = Ax [p] ;
00047 if (i > k && w [i] < k)
00048 {
00049 Vi [vnz++] = i ;
00050 w [i] = k ;
00051 }
00052 }
00053 for (p = top ; p < n ; p++)
00054 {
00055 i = s [p] ;
00056 cs_happly (V, i, Beta [i], x) ;
00057 Ri [rnz] = i ;
00058 Rx [rnz++] = x [i] ;
00059 x [i] = 0 ;
00060 if (parent [i] == k) vnz = cs_scatter (V, i, 0, w, NULL, k, V, vnz);
00061 }
00062 for (p = p1 ; p < vnz ; p++)
00063 {
00064 Vx [p] = x [Vi [p]] ;
00065 x [Vi [p]] = 0 ;
00066 }
00067 Ri [rnz] = k ;
00068 Rx [rnz++] = cs_house (Vx+p1, Beta+k, vnz-p1) ;
00069 }
00070 Rp [n] = rnz ;
00071 Vp [n] = vnz ;
00072 return (cs_ndone (N, NULL, w, x, 1)) ;
00073 }