Go to the documentation of this file.00001 #include "cs.h"
00002
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 = cs_malloc (m+n, sizeof (int)) ;
00008 S->leftmost = leftmost = cs_malloc (m, sizeof (int)) ;
00009 w = cs_malloc (m+3*n, sizeof (int)) ;
00010 if (!pinv || !w || !leftmost)
00011 {
00012 cs_free (w) ;
00013 return (0) ;
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 ;
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 ;
00025 }
00026 }
00027 for (i = m-1 ; i >= 0 ; i--)
00028 {
00029 pinv [i] = -1 ;
00030 k = leftmost [i] ;
00031 if (k == -1) continue ;
00032 if (nque [k]++ == 0) tail [k] = i ;
00033 next [i] = head [k] ;
00034 head [k] = i ;
00035 }
00036 S->lnz = 0 ;
00037 S->m2 = m ;
00038 for (k = 0 ; k < n ; k++)
00039 {
00040 i = head [k] ;
00041 S->lnz++ ;
00042 if (i < 0) i = S->m2++ ;
00043 pinv [i] = k ;
00044 if (--nque [k] <= 0) continue ;
00045 S->lnz += nque [k] ;
00046 if ((pa = parent [k]) != -1)
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
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) ;
00065 n = A->n ;
00066 S = cs_calloc (1, sizeof (css)) ;
00067 if (!S) return (NULL) ;
00068 S->q = cs_amd (order, A) ;
00069 if (order && !S->q) return (cs_sfree (S)) ;
00070 if (qr)
00071 {
00072 cs *C = order ? cs_permute (A, NULL, S->q, 0) : ((cs *) A) ;
00073 S->parent = cs_etree (C, 1) ;
00074 post = cs_post (S->parent, n) ;
00075 S->cp = cs_counts (C, S->parent, post, 1) ;
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 ;
00080 if (order) cs_spfree (C) ;
00081 }
00082 else
00083 {
00084 S->unz = 4*(A->p [n]) + n ;
00085 S->lnz = S->unz ;
00086 }
00087 return (ok ? S : cs_sfree (S)) ;
00088 }