cs_qr.c
Go to the documentation of this file.
00001 #include "cs.h"
00002 /* sparse QR factorization [V,beta,pinv,R] = qr (A) */
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)) ;            /* get int workspace */
00015     x = cs_malloc (m2, sizeof (double)) ;           /* get double workspace */
00016     N = cs_calloc (1, sizeof (csn)) ;               /* allocate result */
00017     if (!w || !x || !N) return (cs_ndone (N, NULL, w, x, 0)) ;
00018     s = w + m2 ;                                    /* s is size n */
00019     for (k = 0 ; k < m2 ; k++) x [k] = 0 ;          /* clear workspace x */
00020     N->L = V = cs_spalloc (m2, n, vnz, 1, 0) ;      /* allocate result V */
00021     N->U = R = cs_spalloc (m2, n, rnz, 1, 0) ;      /* allocate result R */
00022     N->B = Beta = cs_malloc (n, sizeof (double)) ;  /* allocate result Beta */
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 ; /* clear w, to mark nodes */
00027     rnz = 0 ; vnz = 0 ;
00028     for (k = 0 ; k < n ; k++)               /* compute V and R */
00029     {
00030         Rp [k] = rnz ;                      /* R(:,k) starts here */
00031         Vp [k] = p1 = vnz ;                 /* V(:,k) starts here */
00032         w [k] = k ;                         /* add V(k,k) to pattern of V */
00033         Vi [vnz++] = k ;
00034         top = n ;
00035         col = q ? q [k] : k ;
00036         for (p = Ap [col] ; p < Ap [col+1] ; p++)   /* find R(:,k) pattern */
00037         {
00038             i = leftmost [Ai [p]] ;         /* i = min(find(A(i,q))) */
00039             for (len = 0 ; w [i] != k ; i = parent [i]) /* traverse up to k */
00040             {
00041                 s [len++] = i ;
00042                 w [i] = k ;
00043             }
00044             while (len > 0) s [--top] = s [--len] ; /* push path on stack */
00045             i = pinv [Ai [p]] ;             /* i = permuted row of A(:,col) */
00046             x [i] = Ax [p] ;                /* x (i) = A(:,col) */
00047             if (i > k && w [i] < k)         /* pattern of V(:,k) = x (k+1:m) */
00048             {
00049                 Vi [vnz++] = i ;            /* add i to pattern of V(:,k) */
00050                 w [i] = k ;
00051             }
00052         }
00053         for (p = top ; p < n ; p++) /* for each i in pattern of R(:,k) */
00054         {
00055             i = s [p] ;                     /* R(i,k) is nonzero */
00056             cs_happly (V, i, Beta [i], x) ; /* apply (V(i),Beta(i)) to x */
00057             Ri [rnz] = i ;                  /* R(i,k) = x(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++)        /* gather V(:,k) = x */
00063         {
00064             Vx [p] = x [Vi [p]] ;
00065             x [Vi [p]] = 0 ;
00066         }
00067         Ri [rnz] = k ;                     /* R(k,k) = norm (x) */
00068         Rx [rnz++] = cs_house (Vx+p1, Beta+k, vnz-p1) ; /* [v,beta]=house(x) */
00069     }
00070     Rp [n] = rnz ;                          /* finalize R */
00071     Vp [n] = vnz ;                          /* finalize V */
00072     return (cs_ndone (N, NULL, w, x, 1)) ;  /* success */
00073 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines


hogman_minimal
Author(s): Maintained by Juergen Sturm
autogenerated on Wed Dec 26 2012 15:36:47