ufsparse/patch.UFsparse

732 lines
18 KiB
Text

*** UMFPACK/Source/umf_solve.c.orig 2005-02-01 23:53:11.815116040 +0100
--- UMFPACK/Source/umf_solve.c 2005-02-02 00:01:22.904346651 +0100
***************
*** 79,85 ****
double *Z2, *Y, *B2, *Rs ;
Int *Rperm, *Cperm, i, n, p, step, j, nz, status, p2, do_scale ;
#ifdef COMPLEX
! Int split ;
#endif
#ifndef NRECIPROCAL
Int do_recip = Numeric->do_recip ;
--- 79,86 ----
double *Z2, *Y, *B2, *Rs ;
Int *Rperm, *Cperm, i, n, p, step, j, nz, status, p2, do_scale ;
#ifdef COMPLEX
! Int AXsplit ;
! Int Bsplit ;
#endif
#ifndef NRECIPROCAL
Int do_recip = Numeric->do_recip ;
***************
*** 141,147 ****
return (UMFPACK_ERROR_argument_missing) ;
}
/* A, B, and X in split format if Az, Bz, and Xz present */
! split = SPLIT (Az) && SPLIT (Bz) && SPLIT (Xz) ;
Z = (Entry *) (SolveWork + 4*n) ; /* Entry Z [0..n-1] */
S = (Entry *) (SolveWork + 6*n) ; /* Entry S [0..n-1] */
Y = (double *) (SolveWork + 8*n) ; /* double Y [0..n-1] */
--- 142,148 ----
return (UMFPACK_ERROR_argument_missing) ;
}
/* A, B, and X in split format if Az, Bz, and Xz present */
! AXsplit = SPLIT (Az) || SPLIT(Xz);
Z = (Entry *) (SolveWork + 4*n) ; /* Entry Z [0..n-1] */
S = (Entry *) (SolveWork + 6*n) ; /* Entry S [0..n-1] */
Y = (double *) (SolveWork + 8*n) ; /* double Y [0..n-1] */
***************
*** 150,159 ****
}
else
{
! /* A is ignored, only look at X and B for split/packed cases */
! split = SPLIT (Bz) && SPLIT (Xz) ;
}
! if (split)
{
X = (Entry *) (SolveWork + 2*n) ; /* Entry X [0..n-1] */
}
--- 151,162 ----
}
else
{
! /* A is ignored, only look at X for split/packed cases */
! AXsplit = SPLIT(Xz);
}
! Bsplit = SPLIT (Bz);
!
! if (AXsplit)
{
X = (Entry *) (SolveWork + 2*n) ; /* Entry X [0..n-1] */
}
***************
*** 209,215 ****
for (p = 0 ; p < p2 ; p++)
{
/* Y [Ai [p]] += ABS (Ax [p]) ; */
! ASSIGN (aij, Ax, Az, p, split) ;
ABS (d, aij) ;
Y [Ai [p]] += d ;
}
--- 212,218 ----
for (p = 0 ; p < p2 ; p++)
{
/* Y [Ai [p]] += ABS (Ax [p]) ; */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
ABS (d, aij) ;
Y [Ai [p]] += d ;
}
***************
*** 219,225 ****
for (i = 0 ; i < n ; i++)
{
/* B2 [i] = ABS (B [i]) ; */
! ASSIGN (bi, Bx, Bz, i, split) ;
ABS (B2 [i], bi) ;
}
--- 222,228 ----
for (i = 0 ; i < n ; i++)
{
/* B2 [i] = ABS (B [i]) ; */
! ASSIGN (bi, Bx, Bz, i, Bsplit) ;
ABS (B2 [i], bi) ;
}
***************
*** 276,282 ****
/* multiply by the scale factors */
for (i = 0 ; i < n ; i++)
{
! ASSIGN (X [i], Bx, Bz, i, split) ;
SCALE (X [i], Rs [i]) ;
}
}
--- 279,285 ----
/* multiply by the scale factors */
for (i = 0 ; i < n ; i++)
{
! ASSIGN (X [i], Bx, Bz, i, Bsplit) ;
SCALE (X [i], Rs [i]) ;
}
}
***************
*** 286,292 ****
/* divide by the scale factors */
for (i = 0 ; i < n ; i++)
{
! ASSIGN (X [i], Bx, Bz, i, split) ;
SCALE_DIV (X [i], Rs [i]) ;
}
}
--- 289,295 ----
/* divide by the scale factors */
for (i = 0 ; i < n ; i++)
{
! ASSIGN (X [i], Bx, Bz, i, Bsplit) ;
SCALE_DIV (X [i], Rs [i]) ;
}
}
***************
*** 302,308 ****
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [Rperm [i]] ; */
! ASSIGN (W [i], Bx, Bz, Rperm [i], split) ;
}
}
}
--- 305,311 ----
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [Rperm [i]] ; */
! ASSIGN (W [i], Bx, Bz, Rperm [i], Bsplit) ;
}
}
}
***************
*** 311,317 ****
for (i = 0 ; i < n ; i++)
{
/* Z [i] = B [i] ; */
! ASSIGN (Z [i], Bx, Bz, i, split) ;
}
flops += MULTSUB_FLOPS * nz ;
for (i = 0 ; i < n ; i++)
--- 314,320 ----
for (i = 0 ; i < n ; i++)
{
/* Z [i] = B [i] ; */
! ASSIGN (Z [i], Bx, Bz, i, Bsplit) ;
}
flops += MULTSUB_FLOPS * nz ;
for (i = 0 ; i < n ; i++)
***************
*** 321,327 ****
for (p = Ap [i] ; p < p2 ; p++)
{
/* Z [Ai [p]] -= Ax [p] * xi ; */
! ASSIGN (aij, Ax, Az, p, split) ;
MULT_SUB (Z [Ai [p]], aij, xi) ;
}
}
--- 324,330 ----
for (p = Ap [i] ; p < p2 ; p++)
{
/* Z [Ai [p]] -= Ax [p] * xi ; */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
MULT_SUB (Z [Ai [p]], aij, xi) ;
}
}
***************
*** 390,396 ****
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [i] ; */
! ASSIGN (W [i], Bx, Bz, i, split) ;
Z2 [i] = 0. ;
}
flops += (MULT_FLOPS + DECREMENT_FLOPS + ABS_FLOPS + 1) * nz ;
--- 393,399 ----
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [i] ; */
! ASSIGN (W [i], Bx, Bz, i, Bsplit) ;
Z2 [i] = 0. ;
}
flops += (MULT_FLOPS + DECREMENT_FLOPS + ABS_FLOPS + 1) * nz ;
***************
*** 403,409 ****
i = Ai [p] ;
/* axx = Ax [p] * xj ; */
! ASSIGN (aij, Ax, Az, p, split) ;
MULT (axx, aij, xj) ;
/* W [i] -= axx ; */
--- 406,412 ----
i = Ai [p] ;
/* axx = Ax [p] * xj ; */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
MULT (axx, aij, xj) ;
/* W [i] -= axx ; */
***************
*** 493,499 ****
/* yi += ABS (Ax [p]) * Rs [Ai [p]] ; */
/* note that abs (aij) is the same as
* abs (conj (aij)) */
! ASSIGN (aij, Ax, Az, p, split) ;
ABS (d, aij) ;
yi += (d * Rs [Ai [p]]) ;
}
--- 496,502 ----
/* yi += ABS (Ax [p]) * Rs [Ai [p]] ; */
/* note that abs (aij) is the same as
* abs (conj (aij)) */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
ABS (d, aij) ;
yi += (d * Rs [Ai [p]]) ;
}
***************
*** 513,519 ****
/* yi += ABS (Ax [p]) / Rs [Ai [p]] ; */
/* note that abs (aij) is the same as
* abs (conj (aij)) */
! ASSIGN (aij, Ax, Az, p, split) ;
ABS (d, aij) ;
yi += (d / Rs [Ai [p]]) ;
}
--- 516,522 ----
/* yi += ABS (Ax [p]) / Rs [Ai [p]] ; */
/* note that abs (aij) is the same as
* abs (conj (aij)) */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
ABS (d, aij) ;
yi += (d / Rs [Ai [p]]) ;
}
***************
*** 534,540 ****
/* yi += ABS (Ax [p]) ; */
/* note that abs (aij) is the same as
* abs (conj (aij)) */
! ASSIGN (aij, Ax, Az, p, split) ;
ABS (d, aij) ;
yi += d ;
}
--- 537,543 ----
/* yi += ABS (Ax [p]) ; */
/* note that abs (aij) is the same as
* abs (conj (aij)) */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
ABS (d, aij) ;
yi += d ;
}
***************
*** 546,552 ****
for (i = 0 ; i < n ; i++)
{
/* B2 [i] = ABS (B [i]) ; */
! ASSIGN (bi, Bx, Bz, i, split) ;
ABS (B2 [i], bi) ;
}
--- 549,555 ----
for (i = 0 ; i < n ; i++)
{
/* B2 [i] = ABS (B [i]) ; */
! ASSIGN (bi, Bx, Bz, i, Bsplit) ;
ABS (B2 [i], bi) ;
}
***************
*** 568,574 ****
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [Cperm [i]] ; */
! ASSIGN (W [i], Bx, Bz, Cperm [i], split) ;
}
}
else
--- 571,577 ----
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [Cperm [i]] ; */
! ASSIGN (W [i], Bx, Bz, Cperm [i], Bsplit) ;
}
}
else
***************
*** 577,583 ****
for (i = 0 ; i < n ; i++)
{
/* Z [i] = B [i] ; */
! ASSIGN (Z [i], Bx, Bz, i, split) ;
}
flops += MULTSUB_FLOPS * nz ;
for (i = 0 ; i < n ; i++)
--- 580,586 ----
for (i = 0 ; i < n ; i++)
{
/* Z [i] = B [i] ; */
! ASSIGN (Z [i], Bx, Bz, i, Bsplit) ;
}
flops += MULTSUB_FLOPS * nz ;
for (i = 0 ; i < n ; i++)
***************
*** 587,593 ****
for (p = Ap [i] ; p < p2 ; p++)
{
/* zi -= conjugate (Ax [p]) * X [Ai [p]] ; */
! ASSIGN (aij, Ax, Az, p, split) ;
MULT_SUB_CONJ (zi, X [Ai [p]], aij) ;
}
Z [i] = zi ;
--- 590,596 ----
for (p = Ap [i] ; p < p2 ; p++)
{
/* zi -= conjugate (Ax [p]) * X [Ai [p]] ; */
! ASSIGN (aij, Ax, Az, p, Bsplit) ;
MULT_SUB_CONJ (zi, X [Ai [p]], aij) ;
}
Z [i] = zi ;
***************
*** 696,708 ****
for (i = 0 ; i < n ; i++)
{
/* wi = B [i] ; */
! ASSIGN (wi, Bx, Bz, i, split) ;
z2i = 0. ;
p2 = Ap [i+1] ;
for (p = Ap [i] ; p < p2 ; p++)
{
/* axx = conjugate (Ax [p]) * X [Ai [p]] ; */
! ASSIGN (aij, Ax, Az, p, split) ;
MULT_CONJ (axx, X [Ai [p]], aij) ;
/* wi -= axx ; */
--- 699,711 ----
for (i = 0 ; i < n ; i++)
{
/* wi = B [i] ; */
! ASSIGN (wi, Bx, Bz, i, Bsplit) ;
z2i = 0. ;
p2 = Ap [i+1] ;
for (p = Ap [i] ; p < p2 ; p++)
{
/* axx = conjugate (Ax [p]) * X [Ai [p]] ; */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
MULT_CONJ (axx, X [Ai [p]], aij) ;
/* wi -= axx ; */
***************
*** 766,772 ****
/* yi += ABS (Ax [p]) * Rs [Ai [p]] ; */
/* note that A.' is the array transpose,
* so no conjugate */
! ASSIGN (aij, Ax, Az, p, split) ;
ABS (d, aij) ;
yi += (d * Rs [Ai [p]]) ;
}
--- 769,775 ----
/* yi += ABS (Ax [p]) * Rs [Ai [p]] ; */
/* note that A.' is the array transpose,
* so no conjugate */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
ABS (d, aij) ;
yi += (d * Rs [Ai [p]]) ;
}
***************
*** 786,792 ****
/* yi += ABS (Ax [p]) / Rs [Ai [p]] ; */
/* note that A.' is the array transpose,
* so no conjugate */
! ASSIGN (aij, Ax, Az, p, split) ;
ABS (d, aij) ;
yi += (d / Rs [Ai [p]]) ;
}
--- 789,795 ----
/* yi += ABS (Ax [p]) / Rs [Ai [p]] ; */
/* note that A.' is the array transpose,
* so no conjugate */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
ABS (d, aij) ;
yi += (d / Rs [Ai [p]]) ;
}
***************
*** 807,813 ****
/* yi += ABS (Ax [p]) */
/* note that A.' is the array transpose,
* so no conjugate */
! ASSIGN (aij, Ax, Az, p, split) ;
ABS (d, aij) ;
yi += d ;
}
--- 810,816 ----
/* yi += ABS (Ax [p]) */
/* note that A.' is the array transpose,
* so no conjugate */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
ABS (d, aij) ;
yi += d ;
}
***************
*** 819,825 ****
for (i = 0 ; i < n ; i++)
{
/* B2 [i] = ABS (B [i]) ; */
! ASSIGN (bi, Bx, Bz, i, split) ;
ABS (B2 [i], bi) ;
}
--- 822,828 ----
for (i = 0 ; i < n ; i++)
{
/* B2 [i] = ABS (B [i]) ; */
! ASSIGN (bi, Bx, Bz, i, Bsplit) ;
ABS (B2 [i], bi) ;
}
***************
*** 841,847 ****
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [Cperm [i]] ; */
! ASSIGN (W [i], Bx, Bz, Cperm [i], split) ;
}
}
else
--- 844,850 ----
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [Cperm [i]] ; */
! ASSIGN (W [i], Bx, Bz, Cperm [i], Bsplit) ;
}
}
else
***************
*** 850,856 ****
for (i = 0 ; i < n ; i++)
{
/* Z [i] = B [i] ; */
! ASSIGN (Z [i], Bx, Bz, i, split) ;
}
flops += MULTSUB_FLOPS * nz ;
for (i = 0 ; i < n ; i++)
--- 853,859 ----
for (i = 0 ; i < n ; i++)
{
/* Z [i] = B [i] ; */
! ASSIGN (Z [i], Bx, Bz, i, Bsplit) ;
}
flops += MULTSUB_FLOPS * nz ;
for (i = 0 ; i < n ; i++)
***************
*** 860,866 ****
for (p = Ap [i] ; p < p2 ; p++)
{
/* zi -= Ax [p] * X [Ai [p]] ; */
! ASSIGN (aij, Ax, Az, p, split) ;
MULT_SUB (zi, aij, X [Ai [p]]) ;
}
Z [i] = zi ;
--- 863,869 ----
for (p = Ap [i] ; p < p2 ; p++)
{
/* zi -= Ax [p] * X [Ai [p]] ; */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
MULT_SUB (zi, aij, X [Ai [p]]) ;
}
Z [i] = zi ;
***************
*** 969,981 ****
for (i = 0 ; i < n ; i++)
{
/* wi = B [i] ; */
! ASSIGN (wi, Bx, Bz, i, split) ;
z2i = 0. ;
p2 = Ap [i+1] ;
for (p = Ap [i] ; p < p2 ; p++)
{
/* axx = Ax [p] * X [Ai [p]] ; */
! ASSIGN (aij, Ax, Az, p, split) ;
MULT (axx, aij, X [Ai [p]]) ;
/* wi -= axx ; */
--- 972,984 ----
for (i = 0 ; i < n ; i++)
{
/* wi = B [i] ; */
! ASSIGN (wi, Bx, Bz, i, Bsplit) ;
z2i = 0. ;
p2 = Ap [i+1] ;
for (p = Ap [i] ; p < p2 ; p++)
{
/* axx = Ax [p] * X [Ai [p]] ; */
! ASSIGN (aij, Ax, Az, p, AXsplit) ;
MULT (axx, aij, X [Ai [p]]) ;
/* wi -= axx ; */
***************
*** 1011,1017 ****
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [Rperm [i]] ; */
! ASSIGN (X [i], Bx, Bz, Rperm [i], split) ;
}
flops = UMF_lsolve (Numeric, X, Pattern) ;
status = UMFPACK_OK ;
--- 1014,1020 ----
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [Rperm [i]] ; */
! ASSIGN (X [i], Bx, Bz, Rperm [i], Bsplit) ;
}
flops = UMF_lsolve (Numeric, X, Pattern) ;
status = UMFPACK_OK ;
***************
*** 1027,1033 ****
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, split) ;
}
flops = UMF_lsolve (Numeric, X, Pattern) ;
status = UMFPACK_OK ;
--- 1030,1036 ----
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, Bsplit) ;
}
flops = UMF_lsolve (Numeric, X, Pattern) ;
status = UMFPACK_OK ;
***************
*** 1043,1049 ****
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [i] ; */
! ASSIGN (W [i], Bx, Bz, i, split) ;
}
flops = UMF_lhsolve (Numeric, W, Pattern) ;
for (i = 0 ; i < n ; i++)
--- 1046,1052 ----
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [i] ; */
! ASSIGN (W [i], Bx, Bz, i, Bsplit) ;
}
flops = UMF_lhsolve (Numeric, W, Pattern) ;
for (i = 0 ; i < n ; i++)
***************
*** 1063,1069 ****
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [i] ; */
! ASSIGN (W [i], Bx, Bz, i, split) ;
}
flops = UMF_ltsolve (Numeric, W, Pattern) ;
for (i = 0 ; i < n ; i++)
--- 1066,1072 ----
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [i] ; */
! ASSIGN (W [i], Bx, Bz, i, Bsplit) ;
}
flops = UMF_ltsolve (Numeric, W, Pattern) ;
for (i = 0 ; i < n ; i++)
***************
*** 1083,1089 ****
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, split) ;
}
flops = UMF_lhsolve (Numeric, X, Pattern) ;
status = UMFPACK_OK ;
--- 1086,1092 ----
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, Bsplit) ;
}
flops = UMF_lhsolve (Numeric, X, Pattern) ;
status = UMFPACK_OK ;
***************
*** 1099,1105 ****
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, split) ;
}
flops = UMF_ltsolve (Numeric, X, Pattern) ;
status = UMFPACK_OK ;
--- 1102,1108 ----
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, Bsplit) ;
}
flops = UMF_ltsolve (Numeric, X, Pattern) ;
status = UMFPACK_OK ;
***************
*** 1115,1121 ****
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [i] ; */
! ASSIGN (W [i], Bx, Bz, i, split) ;
}
flops = UMF_usolve (Numeric, W, Pattern) ;
for (i = 0 ; i < n ; i++)
--- 1118,1124 ----
for (i = 0 ; i < n ; i++)
{
/* W [i] = B [i] ; */
! ASSIGN (W [i], Bx, Bz, i, Bsplit) ;
}
flops = UMF_usolve (Numeric, W, Pattern) ;
for (i = 0 ; i < n ; i++)
***************
*** 1134,1140 ****
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, split) ;
}
flops = UMF_usolve (Numeric, X, Pattern) ;
--- 1137,1143 ----
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, Bsplit) ;
}
flops = UMF_usolve (Numeric, X, Pattern) ;
***************
*** 1149,1155 ****
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [Cperm [i]] ; */
! ASSIGN (X [i], Bx, Bz, Cperm [i], split) ;
}
flops = UMF_uhsolve (Numeric, X, Pattern) ;
--- 1152,1158 ----
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [Cperm [i]] ; */
! ASSIGN (X [i], Bx, Bz, Cperm [i], Bsplit) ;
}
flops = UMF_uhsolve (Numeric, X, Pattern) ;
***************
*** 1164,1170 ****
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [Cperm [i]] ; */
! ASSIGN (X [i], Bx, Bz, Cperm [i], split) ;
}
flops = UMF_utsolve (Numeric, X, Pattern) ;
--- 1167,1173 ----
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [Cperm [i]] ; */
! ASSIGN (X [i], Bx, Bz, Cperm [i], Bsplit) ;
}
flops = UMF_utsolve (Numeric, X, Pattern) ;
***************
*** 1179,1185 ****
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, split) ;
}
flops = UMF_uhsolve (Numeric, X, Pattern) ;
--- 1182,1188 ----
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, Bsplit) ;
}
flops = UMF_uhsolve (Numeric, X, Pattern) ;
***************
*** 1194,1200 ****
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, split) ;
}
flops = UMF_utsolve (Numeric, X, Pattern) ;
--- 1197,1203 ----
for (i = 0 ; i < n ; i++)
{
/* X [i] = B [i] ; */
! ASSIGN (X [i], Bx, Bz, i, Bsplit) ;
}
flops = UMF_utsolve (Numeric, X, Pattern) ;
***************
*** 1206,1212 ****
#ifdef COMPLEX
/* copy the solution back, from Entry X [ ] to double Xx [ ] and Xz [ ] */
! if (split)
{
for (i = 0 ; i < n ; i++)
{
--- 1209,1215 ----
#ifdef COMPLEX
/* copy the solution back, from Entry X [ ] to double Xx [ ] and Xz [ ] */
! if (AXsplit)
{
for (i = 0 ; i < n ; i++)
{