201 int p, j, nz = 0, anz, *Cp, *Cj, *Bp, m, n, bnz, *w, values ;
202 double *x, *Bx, *Cx ;
205 if ( (A->
m != B->
m) || (A->
n != B->
n) )
return (NULL);
206 m = A->
m ; anz = A->
p [m] ;
207 n = B->
n ; Bp = B->
p ; Bx = B->
x ; bnz = Bp [m] ;
208 w = (
int*)calloc (
CS_MAX(n,1),
sizeof (int)) ;
209 values = (A->
x != NULL) && (Bx != NULL) ;
210 x = values ? (
double*)malloc (n *
sizeof (
double)) : NULL ;
212 if (!C || !w || (values && !x))
return (
csr_done (C, w, x, 0)) ;
213 Cp = C->
p ; Cj = C->
j ; Cx = C->
x ;
214 for (j = 0 ; j < n ; j++)
219 if (values)
for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Cj [p]] ;
254 int p, j, nz = 0, anz, *Cp, *Cj, *Ap, m, k, n, bnz, *w, values, *Aj;
255 double *x, *Ax, *Cx ;
259 if (k != B->
m )
return(NULL);
260 m = A->
m ; anz = A->
p [A->
m] ;
261 n = B->
n ; Ap = A->
p ; Aj = A->
j ; Ax = A->
x ; bnz = B->
p [k] ;
262 w = (
int*)calloc (
CS_MAX(n,1),
sizeof (int)) ;
263 values = (Ax != NULL) && (B->
x != NULL) ;
264 x = values ? (
double*)malloc (n *
sizeof (
double)) : NULL ;
266 if (!C || !w || (values && !x))
return (
csr_done (C, w, x, 0)) ;
268 for (j = 0 ; j < m ; j++)
274 Cj = C->
j ; Cx = C->
x ;
276 for (p = Ap [j] ; p < Ap [j+1] ; p++)
278 nz =
csr_scatter (B, Aj [p], Ax ? Ax [p] : 1, w, x, j+1, C, nz) ;
280 if (values)
for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Cj [p]] ;
407 double pivot, *Lx, *Ux, *x, a, t;
408 int *Lp, *Lj, *Up, *Uj, *pinv, *xi, *q, n, ipiv, k, top, p, i, row, lnz, unz;
411 if (!
CS_CSC (A) ) printf(
" error csrlu: A not csc\n");
412 if (!
CS_CSC (A) || !S)
return (NULL) ;
414 if (n != A->
m)
return (NULL) ;
415 q = S->
q ; lnz = (int)S->
lnz ; unz = (int)S->
unz ;
416 x = (
double*)malloc(n *
sizeof (
double)) ;
417 xi = (
int*)malloc (2*n *
sizeof (
int)) ;
419 if (!(A) ) printf(
" error csrlu: allocation of N failed\n");
420 if (!x || !xi || !
N)
return (
csr_ndone (
N, NULL, xi, x, 0)) ;
424 N->pinv = pinv = (
int*)malloc (n *
sizeof (
int)) ;
425 N->perm = (
int*)malloc (n *
sizeof (
int)) ;
426 if (!L || !U || !pinv)
return (
csr_ndone (
N, NULL, xi, x, 0)) ;
427 Lp = L->
p ; Up = U->
p ;
428 for (i = 0 ; i < n ; i++) x [i] = 0 ;
429 for (i = 0 ; i < n ; i++) pinv [i] = -1 ;
430 for (k = 1 ; k <= n ; k++) Up [k] = 0 ;
431 for (k = 1 ; k <= n ; k++) Lp [k] = 0 ;
437 for (k = 0 ; k < n ; k++)
447 Lj = L->
j ; Lx = L->
x ; Uj = U->
j ; Ux = U->
x ;
448 row = q ? (q [k]) : k ;
451 printf(
"--------------------------------\n");
452 printf(
" %d spsolve row=%d \n",k, row);
453 printf(
" pinv = %d %d %d %d \n", pinv[0], pinv[1], pinv[2], pinv[3]);
456 if( debug > 1 ) printf(
"top=%d x = %g %g %g %g \n", top,x[0],x[1],x[2],x[3]);
460 for (p = top ; p < n ; p++)
465 if ((t = fabs (x [i])) > a)
473 Lj [lnz] = pinv [i] ;
477 if (ipiv == -1 || a <= 0)
return (
csr_ndone (
N, NULL, xi, x, 0)) ;
478 if (pinv [row] < 0 && fabs (x [row]) >= a*
tol) ipiv = row;
482 if( debug > 1 ) { printf (
"L:") ;
csr_print (L, 0) ; }
488 for (p = top ; p < n ; p++)
494 Ux [unz++] = x [i] / pivot ;
501 printf(
"------------------------------------\n");
506 if( debug ) { printf (
"L:") ;
csr_print (L, 0) ; }
509 for (p = 0 ; p < unz ; p++) Uj [p] = pinv [Uj [p]] ;
513 if( debug ) { printf (
"U:") ;
csr_print (U, 0) ; }
524 int i, I, p, q, px, top, n, *Gp, *Gj, *Bp, *Bj ;
527 if (!
CS_CSC (G) || !
CS_CSC (B) || !xi || !x)
return (-1) ;
528 Gp = G->
p ; Gj = G->
j ; Gx = G->
x ; n = G->
n ;
529 Bp = B->
p ; Bj = B->
j ; Bx = B->
x ;
531 for (p = top ; p < n ; p++) x [xi [p]] = 0 ;
532 for (p = Bp [k] ; p < Bp [k+1] ; p++) x [Bj [p]] = Bx [p] ;
534 printf(
"solve k=%d x= %g %g %g %g \n", k, x[0], x[1], x[2], x[3]);
536 printf(
"top=%d xi= %d %d %d %d \n", top , xi[0], xi[1], xi[2], xi[3]);
537 for (px = top ; px < n ; px++)
540 I = pinv ? (pinv [i]) : i ;
541 if (I < 0) continue ;
543 x [i] /= Gx [up ? (Gp [I]) : (Gp [I+1]-1)] ;
544 p = up ? (Gp [I]+1) : (Gp [I]) ;
545 q = up ? (Gp [I+1]) : (Gp [I+1]-1) ;
549 printf(
"%d %d solve %d %g %g \n", px, i ,p, Gx [p] , x [Gj [p]] );
551 x [Gj[p]] -= Gx [p] * x [i] ;
554 printf(
" x= %g %g %g %g \n", x[0], x[1], x[2], x[3]);
586 int *Cp, *Cj, *last, *W, *len, *nv, *next, *P, *head, *elen, *degree, *w,
587 *hhead, *ATp, *ATj, d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
588 k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
589 ok, cnz, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, n, m, t ;
592 if (!
CS_CSC (A) || order <= 0 || order > 3)
return (NULL) ;
594 if (!AT)
return (NULL) ;
595 m = A->
m ; n = A->
n ;
596 if ( n != m)
return(NULL);
597 dense = (int)
CS_MAX (16, 10 * sqrt ((
double) n)) ;
598 dense =
CS_MIN (n-2, dense) ;
599 if (order == 1 && n == m)
607 for (p2 = 0, j = 0 ; j < m ; j++)
611 if (ATp [j+1] - p > dense) continue ;
612 for ( ; p < ATp [j+1] ; p++) ATj [p2++] = ATj [p] ;
624 if (!C)
return (NULL) ;
628 P = (
int*)malloc (
CS_MAX(n+1,1) *
sizeof (int)) ;
629 W = (
int*)malloc (
CS_MAX(8*(n+1),1) *
sizeof (
int)) ;
630 t = cnz + cnz/5 + 2*n ;
632 len = W ; nv = W + (n+1) ; next = W + 2*(n+1) ;
633 head = W + 3*(n+1) ; elen = W + 4*(n+1) ; degree = W + 5*(n+1) ;
634 w = W + 6*(n+1) ; hhead = W + 7*(n+1) ;
637 for (k = 0 ; k < n ; k++) len [k] = Cp [k+1] - Cp [k] ;
641 for (i = 0 ; i <= n ; i++)
650 degree [i] = len [i] ;
657 for (i = 0 ; i < n ; i++)
677 if (head [d] != -1) last [head [d]] = i ;
678 next [i] = head [d] ;
685 for (k = -1 ; mindeg < n && (k = head [mindeg]) == -1 ; mindeg++) ;
686 if (next [k] != -1) last [next [k]] = -1 ;
687 head [mindeg] = next [k] ;
692 if (elenk > 0 && cnz + mindeg >= nzmax)
694 for (j = 0 ; j < n ; j++)
696 if ((p = Cp [j]) >= 0)
702 for (q = 0, p = 0 ; p < cnz ; )
704 if ((j =
CS_FLIP (Cj [p++])) >= 0)
708 for (k3 = 0 ; k3 < len [j]-1 ; k3++) Cj [q++] = Cj [p++] ;
717 pk1 = (elenk == 0) ? p : cnz ;
719 for (k1 = 1 ; k1 <= elenk + 1 ; k1++)
725 ln = len [k] - elenk ;
733 for (k2 = 1 ; k2 <= ln ; k2++)
736 if ((nvi = nv [i]) <= 0) continue ;
740 if (next [i] != -1) last [next [i]] = last [i] ;
743 next [last [i]] = next [i] ;
747 head [degree [i]] = next [i] ;
756 if (elenk != 0) cnz = pk2 ;
759 len [k] = pk2 - pk1 ;
763 for (pk = pk1 ; pk < pk2 ; pk++)
766 if ((eln = elen [i]) <= 0) continue ;
769 for (p = Cp [i] ; p <= Cp [i] + eln - 1 ; p++)
778 w [e] = degree [e] + wnvi ;
783 for (pk = pk1 ; pk < pk2 ; pk++)
787 p2 = p1 + elen [i] - 1 ;
789 for (h = 0, d = 0, p = p1 ; p <= p2 ; p++)
794 dext = w [e] - mark ;
808 elen [i] = pn - p1 + 1 ;
811 for (p = p2 + 1 ; p < p4 ; p++)
814 if ((nvj = nv [j]) <= 0) continue ;
831 degree [i] =
CS_MIN (degree [i], d) ;
836 len [i] = pn - p1 + 1 ;
838 next [i] = hhead [h] ;
844 lemax =
CS_MAX (lemax, dk) ;
847 for (pk = pk1 ; pk < pk2 ; pk++)
850 if (nv [i] >= 0) continue ;
854 for ( ; i != -1 && next [i] != -1 ; i = next [i], mark++)
858 for (p = Cp [i]+1 ; p <= Cp [i] + ln-1 ; p++) w [Cj [p]] = mark;
860 for (j = next [i] ; j != -1 ; )
862 ok = (len [j] == ln) && (elen [j] == eln) ;
863 for (p = Cp [j] + 1 ; ok && p <= Cp [j] + ln - 1 ; p++)
865 if (w [Cj [p]] != mark) ok = 0 ;
885 for (p = pk1, pk = pk1 ; pk < pk2 ; pk++)
888 if ((nvi = -nv [i]) <= 0) continue ;
890 d = degree [i] + dk - nvi ;
891 d =
CS_MIN (d, n - nel - nvi) ;
892 if (head [d] != -1) last [head [d]] = i ;
893 next [i] = head [d] ;
896 mindeg =
CS_MIN (mindeg, d) ;
901 if ((len [k] = p-pk1) == 0)
906 if (elenk != 0) cnz = p ;
909 for (i = 0 ; i < n ; i++) Cp [i] =
CS_FLIP (Cp [i]) ;
910 for (j = 0 ; j <= n ; j++) head [j] = -1 ;
911 for (j = n ; j >= 0 ; j--)
913 if (nv [j] > 0) continue ;
914 next [j] = head [Cp [j]] ;
917 for (e = n ; e >= 0 ; e--)
919 if (nv [e] <= 0) continue ;
922 next [e] = head [Cp [e]] ;
926 for (k = 0, i = 0 ; i <= n ; i++)
928 if (Cp [i] == -1) k =
csr_tdfs (i, k, head, next, P, w) ;
941 int p, j, m, n, nzmax, nz, nnz, *Ap, *Aj ;
943 if (!A) { printf (
"(null)\n") ;
return (0) ; }
944 m = A->
m ; n = A->
n ; Ap = A->
p ; Aj = A->
j ; Ax = A->
x ;
945 nzmax = A->
nzmax ; nz = A->
nz ; nnz = Ap [m];
950 while ((Ap[m] == 0) && (m > 0))
958 printf (
"%d-by-%d, nzmax: %d nnz: %d, mxnorm: %g\n", m, n, nzmax,
960 for (j = 0 ; j < m ; j++)
962 printf (
" row %d : locations %d to %d\n", j, Ap [j], Ap [j+1]-1);
963 for (p = Ap [j] ; p < Ap [j+1] ; p++)
965 printf (
" %d : %g\n", Aj [p], Ax ? Ax [p] : 1) ;
966 if (brief && p > 20) { printf (
" ...\n") ;
return (1) ; }
970 printf (
"%d-by-%d, zero matrix with nzmax: %d\n", m, n, nzmax);
975 printf (
"triplet: %d-by-%d, nzmax: %d nnz: %d\n", m, n, nzmax, nz) ;
976 for (p = 0 ; p < nz ; p++)
978 printf (
" %d %d : %g\n", Aj [p], Ap [p], Ax ? Ax [p] : 1) ;
979 if (brief && p > 20) { printf (
" ...\n") ;
return (1) ; }