cs_sqr.c
Go to the documentation of this file.
00001 #include "cs.h"
00002 /* compute nnz(V) = S->lnz, S->pinv, S->leftmost, S->m2 from A and S->parent */
00003 static int cs_vcount (const cs *A, css *S)
00004 {
00005     int i, k, p, pa, n = A->n, m = A->m, *Ap = A->p, *Ai = A->i, *next, *head,
00006         *tail, *nque, *pinv, *leftmost, *w, *parent = S->parent ;
00007     S->pinv = pinv = (int*)cs_malloc (m+n, sizeof (int)) ;        /* allocate pinv, */
00008     S->leftmost = leftmost = (int*)cs_malloc (m, sizeof (int)) ;  /* and leftmost */
00009     w = (int*)cs_malloc (m+3*n, sizeof (int)) ;   /* get workspace */
00010     if (!pinv || !w || !leftmost)
00011     {
00012         cs_free (w) ;                       /* pinv and leftmost freed later */
00013         return (0) ;                        /* out of memory */
00014     }
00015     next = w ; head = w + m ; tail = w + m + n ; nque = w + m + 2*n ;
00016     for (k = 0 ; k < n ; k++) head [k] = -1 ;   /* queue k is empty */
00017     for (k = 0 ; k < n ; k++) tail [k] = -1 ;
00018     for (k = 0 ; k < n ; k++) nque [k] = 0 ;
00019     for (i = 0 ; i < m ; i++) leftmost [i] = -1 ;
00020     for (k = n-1 ; k >= 0 ; k--)
00021     {
00022         for (p = Ap [k] ; p < Ap [k+1] ; p++)
00023         {
00024             leftmost [Ai [p]] = k ;         /* leftmost[i] = min(find(A(i,:)))*/
00025         }
00026     }
00027     for (i = m-1 ; i >= 0 ; i--)            /* scan rows in reverse order */
00028     {
00029         pinv [i] = -1 ;                     /* row i is not yet ordered */
00030         k = leftmost [i] ;
00031         if (k == -1) continue ;             /* row i is empty */
00032         if (nque [k]++ == 0) tail [k] = i ; /* first row in queue k */
00033         next [i] = head [k] ;               /* put i at head of queue k */
00034         head [k] = i ;
00035     }
00036     S->lnz = 0 ;
00037     S->m2 = m ;
00038     for (k = 0 ; k < n ; k++)               /* find row permutation and nnz(V)*/
00039     {
00040         i = head [k] ;                      /* remove row i from queue k */
00041         S->lnz++ ;                          /* count V(k,k) as nonzero */
00042         if (i < 0) i = S->m2++ ;            /* add a fictitious row */
00043         pinv [i] = k ;                      /* associate row i with V(:,k) */
00044         if (--nque [k] <= 0) continue ;     /* skip if V(k+1:m,k) is empty */
00045         S->lnz += nque [k] ;                /* nque [k] is nnz (V(k+1:m,k)) */
00046         if ((pa = parent [k]) != -1)        /* move all rows to parent of k */
00047         {
00048             if (nque [pa] == 0) tail [pa] = tail [k] ;
00049             next [tail [k]] = head [pa] ;
00050             head [pa] = next [i] ;
00051             nque [pa] += nque [k] ;
00052         }
00053     }
00054     for (i = 0 ; i < m ; i++) if (pinv [i] < 0) pinv [i] = k++ ;
00055     cs_free (w) ;
00056     return (1) ;
00057 }
00058 
00059 /* symbolic ordering and analysis for QR or LU */
00060 css *cs_sqr (int order, const cs *A, int qr)
00061 {
00062     int n, k, ok = 1, *post ;
00063     css *S ;
00064     if (!CS_CSC (A)) return (NULL) ;        /* check inputs */
00065     n = A->n ;
00066     S = (css*)cs_calloc (1, sizeof (css)) ;       /* allocate result S */
00067     if (!S) return (NULL) ;                 /* out of memory */
00068     S->q = cs_amd (order, A) ;              /* fill-reducing ordering */
00069     if (order && !S->q) return (cs_sfree (S)) ;
00070     if (qr)                                 /* QR symbolic analysis */
00071     {
00072         cs *C = order ? cs_permute (A, NULL, S->q, 0) : ((cs *) A) ;
00073         S->parent = cs_etree (C, 1) ;       /* etree of C'*C, where C=A(:,q) */
00074         post = cs_post (S->parent, n) ;
00075         S->cp = cs_counts (C, S->parent, post, 1) ;  /* col counts chol(C'*C) */
00076         cs_free (post) ;
00077         ok = C && S->parent && S->cp && cs_vcount (C, S) ;
00078         if (ok) for (S->unz = 0, k = 0 ; k < n ; k++) S->unz += S->cp [k] ;
00079         ok = ok && S->lnz >= 0 && S->unz >= 0 ;     /* int overflow guard */
00080         if (order) cs_spfree (C) ;
00081     }
00082     else
00083     {
00084         S->unz = 4*(A->p [n]) + n ;         /* for LU factorization only, */
00085         S->lnz = S->unz ;                   /* guess nnz(L) and nnz(U) */
00086     }
00087     return (ok ? S : cs_sfree (S)) ;        /* return result S */
00088 }


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:58:04