1
2 /* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 File lusol6a
4 lu6sol lu6L lu6Lt lu6U Lu6Ut lu6LD lu6chk
5 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6 26 Apr 2002: lu6 routines put into a separate file.
7 15 Dec 2002: lu6sol modularized via lu6L, lu6Lt, lu6U, lu6Ut.
8 lu6LD implemented to allow solves with LDL' or L|D|L'.
9 15 Dec 2002: Current version of lusol6a.f.
10 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
11
12 /* ==================================================================
13 lu6chk looks at the LU factorization A = L*U.
14 If mode = 1, lu6chk is being called by lu1fac.
15 (Other modes not yet implemented.)
16 ------------------------------------------------------------------
17 The important input parameters are
18
19 lprint = luparm(2)
20 keepLU = luparm(8)
21 Utol1 = parmlu(4)
22 Utol2 = parmlu(5)
23
24 and the significant output parameters are
25
26 inform = luparm(10)
27 nsing = luparm(11)
28 jsing = luparm(12)
29 jumin = luparm(19)
30 Lmax = parmlu(11)
31 Umax = parmlu(12)
32 DUmax = parmlu(13)
33 DUmin = parmlu(14)
34 and w(*).
35
36 Lmax and Umax return the largest elements in L and U.
37 DUmax and DUmin return the largest and smallest diagonals of U
38 (excluding diagonals that are exactly zero).
39 In general, w(j) is set to the maximum absolute element in
40 the j-th column of U. However, if the corresponding diagonal
41 of U is small in absolute terms or relative to w(j)
42 (as judged by the parameters Utol1, Utol2 respectively),
43 then w(j) is changed to - w(j).
44 Thus, if w(j) is not positive, the j-th column of A
45 appears to be dependent on the other columns of A.
46 The number of such columns, and the position of the last one,
47 are returned as nsing and jsing.
48 Note that nrank is assumed to be set already, and is not altered.
49 Typically, nsing will satisfy nrank + nsing = n, but if
50 Utol1 and Utol2 are rather large, nsing > n - nrank may occur.
51 If keepLU = 0,
52 Lmax and Umax are already set by lu1fac.
53 The diagonals of U are in the top of A.
54 Only Utol1 is used in the singularity test to set w(*).
55 inform = 0 if A appears to have full column rank (nsing = 0).
56 inform = 1 otherwise (nsing .gt. 0).
57 ------------------------------------------------------------------
58 00 Jul 1987: Early version.
59 09 May 1988: f77 version.
60 11 Mar 2001: Allow for keepLU = 0.
61 17 Nov 2001: Briefer output for singular factors.
62 05 May 2002: Comma needed in format 1100 (via Kenneth Holmstrom).
63 06 May 2002: With keepLU = 0, diags of U are in natural order.
64 They were not being extracted correctly.
65 23 Apr 2004: TRP can judge singularity better by comparing
66 all diagonals to DUmax.
67 27 Jun 2004: (PEG) Allow write only if nout .gt. 0.
68 ================================================================== */
69 #ifdef UseOld_LU6CHK_20040510
LU6CHK(LUSOLrec * LUSOL,int MODE,int LENA2,int * INFORM)70 void LU6CHK(LUSOLrec *LUSOL, int MODE, int LENA2, int *INFORM)
71 {
72 MYBOOL KEEPLU;
73 int I, J, JUMIN, K, L, L1, L2, LENL, LPRINT, NDEFIC, NRANK;
74 REAL AIJ, DIAG, DUMAX, DUMIN, LMAX, UMAX, UTOL1, UTOL2;
75
76 LPRINT = LUSOL->luparm[LUSOL_IP_PRINTLEVEL];
77 KEEPLU = (MYBOOL) (LUSOL->luparm[LUSOL_IP_KEEPLU]!=0);
78 NRANK = LUSOL->luparm[LUSOL_IP_RANK_U];
79 LENL = LUSOL->luparm[LUSOL_IP_NONZEROS_L];
80 UTOL1 = LUSOL->parmlu[LUSOL_RP_SMALLDIAG_U];
81 UTOL2 = LUSOL->parmlu[LUSOL_RP_EPSDIAG_U];
82 *INFORM = LUSOL_INFORM_LUSUCCESS;
83 LMAX = ZERO;
84 UMAX = ZERO;
85 LUSOL->luparm[LUSOL_IP_SINGULARITIES] = 0;
86 LUSOL->luparm[LUSOL_IP_SINGULARINDEX] = 0;
87 JUMIN = 0;
88 DUMAX = ZERO;
89 DUMIN = LUSOL_BIGNUM;
90
91 #ifdef LUSOLFastClear
92 MEMCLEAR(LUSOL->w+1, LUSOL->n);
93 #else
94 for(I = 1; I <= LUSOL->n; I++)
95 LUSOL->w[I] = ZERO;
96 #endif
97
98 if(KEEPLU) {
99 /* --------------------------------------------------------------
100 Find Lmax.
101 -------------------------------------------------------------- */
102 for(L = (LENA2+1)-LENL; L <= LENA2; L++) {
103 SETMAX(LMAX,fabs(LUSOL->a[L]));
104 }
105 /* --------------------------------------------------------------
106 Find Umax and set w(j) = maximum element in j-th column of U.
107 -------------------------------------------------------------- */
108 for(K = 1; K <= NRANK; K++) {
109 I = LUSOL->ip[K];
110 L1 = LUSOL->locr[I];
111 L2 = (L1+LUSOL->lenr[I])-1;
112 for(L = L1; L <= L2; L++) {
113 J = LUSOL->indr[L];
114 AIJ = fabs(LUSOL->a[L]);
115 SETMAX(LUSOL->w[J],AIJ);
116 SETMAX(UMAX,AIJ);
117 }
118 }
119 /* --------------------------------------------------------------
120 Negate w(j) if the corresponding diagonal of U is
121 too small in absolute terms or relative to the other elements
122 in the same column of U.
123 Also find DUmax and DUmin, the extreme diagonals of U.
124 -------------------------------------------------------------- */
125 for(K = 1; K <= LUSOL->n; K++) {
126 J = LUSOL->iq[K];
127 if(K>NRANK)
128 DIAG = ZERO;
129 else {
130 I = LUSOL->ip[K];
131 L1 = LUSOL->locr[I];
132 DIAG = fabs(LUSOL->a[L1]);
133 SETMAX(DUMAX,DIAG);
134 if(DUMIN>DIAG) {
135 DUMIN = DIAG;
136 JUMIN = J;
137 }
138 }
139 if((DIAG<=UTOL1) || (DIAG<=UTOL2*LUSOL->w[J])) {
140 LUSOL_addSingularity(LUSOL, J, INFORM);
141 LUSOL->w[J] = -LUSOL->w[J];
142 }
143 }
144 LUSOL->parmlu[LUSOL_RP_MAXMULT_L] = LMAX;
145 LUSOL->parmlu[LUSOL_RP_MAXELEM_U] = UMAX;
146 }
147 else {
148 /* --------------------------------------------------------------
149 keepLU = 0.
150 Only diag(U) is stored. Set w(*) accordingly.
151 -------------------------------------------------------------- */
152 for(K = 1; K <= LUSOL->n; K++) {
153 J = LUSOL->iq[K];
154 if(K>NRANK)
155 DIAG = ZERO;
156 else {
157 /* ! diag = abs( diagU(k) ) ! 06 May 2002: Diags are in natural order */
158 DIAG = fabs(LUSOL->diagU[J]);
159 LUSOL->w[J] = DIAG;
160 SETMAX(DUMAX,DIAG);
161 if(DUMIN>DIAG) {
162 DUMIN = DIAG;
163 JUMIN = J;
164 }
165 }
166 if(DIAG<=UTOL1) {
167 LUSOL_addSingularity(LUSOL, J, INFORM);
168 LUSOL->w[J] = -LUSOL->w[J];
169 }
170 }
171 }
172 /* -----------------------------------------------------------------
173 Set output parameters.
174 ----------------------------------------------------------------- */
175 if(JUMIN==0)
176 DUMIN = ZERO;
177 LUSOL->luparm[LUSOL_IP_COLINDEX_DUMIN] = JUMIN;
178 LUSOL->parmlu[LUSOL_RP_MAXELEM_DIAGU] = DUMAX;
179 LUSOL->parmlu[LUSOL_RP_MINELEM_DIAGU] = DUMIN;
180 /* The matrix has been judged singular. */
181 if(LUSOL->luparm[LUSOL_IP_SINGULARITIES]>0) {
182 *INFORM = LUSOL_INFORM_LUSINGULAR;
183 NDEFIC = LUSOL->n-NRANK;
184 if(LPRINT>=LUSOL_MSG_SINGULARITY) {
185 LUSOL_report(LUSOL, 0, "Singular(m%cn) rank:%9d n-rank:%8d nsing:%9d\n",
186 relationChar(LUSOL->m, LUSOL->n),NRANK,NDEFIC,
187 LUSOL->luparm[LUSOL_IP_SINGULARITIES]);
188 }
189 }
190 /* Exit. */
191 LUSOL->luparm[LUSOL_IP_INFORM] = *INFORM;
192 }
193 #else
LU6CHK(LUSOLrec * LUSOL,int MODE,int LENA2,int * INFORM)194 void LU6CHK(LUSOLrec *LUSOL, int MODE, int LENA2, int *INFORM)
195 {
196 MYBOOL KEEPLU, TRP;
197 int I, J, JUMIN, K, L, L1, L2, LENL, LDIAGU, LPRINT, NDEFIC, NRANK;
198 REAL AIJ, DIAG, DUMAX, DUMIN, LMAX, UMAX, UTOL1, UTOL2;
199
200 LPRINT = LUSOL->luparm[LUSOL_IP_PRINTLEVEL];
201 KEEPLU = (MYBOOL) (LUSOL->luparm[LUSOL_IP_KEEPLU] != 0);
202 TRP = (MYBOOL) (LUSOL->luparm[LUSOL_IP_PIVOTTYPE] == LUSOL_PIVMOD_TRP);
203 NRANK = LUSOL->luparm[LUSOL_IP_RANK_U];
204 LENL = LUSOL->luparm[LUSOL_IP_NONZEROS_L];
205 UTOL1 = LUSOL->parmlu[LUSOL_RP_SMALLDIAG_U];
206 UTOL2 = LUSOL->parmlu[LUSOL_RP_EPSDIAG_U];
207 *INFORM = LUSOL_INFORM_LUSUCCESS;
208 LMAX = ZERO;
209 UMAX = ZERO;
210 LUSOL->luparm[LUSOL_IP_SINGULARITIES] = 0;
211 LUSOL->luparm[LUSOL_IP_SINGULARINDEX] = 0;
212 JUMIN = 0;
213 DUMAX = ZERO;
214 DUMIN = LUSOL_BIGNUM;
215
216 #ifdef LUSOLFastClear
217 MEMCLEAR(LUSOL->w+1, LUSOL->n);
218 #else
219 for(I = 1; I <= LUSOL->n; I++)
220 LUSOL->w[I] = ZERO;
221 #endif
222
223 if(KEEPLU) {
224 /* --------------------------------------------------------------
225 Find Lmax.
226 -------------------------------------------------------------- */
227 for(L = (LENA2+1)-LENL; L <= LENA2; L++) {
228 SETMAX(LMAX,fabs(LUSOL->a[L]));
229 }
230 /* --------------------------------------------------------------
231 Find Umax and set w(j) = maximum element in j-th column of U.
232 -------------------------------------------------------------- */
233 for(K = 1; K <= NRANK; K++) {
234 I = LUSOL->ip[K];
235 L1 = LUSOL->locr[I];
236 L2 = (L1+LUSOL->lenr[I])-1;
237 for(L = L1; L <= L2; L++) {
238 J = LUSOL->indr[L];
239 AIJ = fabs(LUSOL->a[L]);
240 SETMAX(LUSOL->w[J],AIJ);
241 SETMAX(UMAX,AIJ);
242 }
243 }
244 LUSOL->parmlu[LUSOL_RP_MAXMULT_L] = LMAX;
245 LUSOL->parmlu[LUSOL_RP_MAXELEM_U] = UMAX;
246 /* --------------------------------------------------------------
247 Find DUmax and DUmin, the extreme diagonals of U.
248 -------------------------------------------------------------- */
249 for(K = 1; K <= NRANK; K++) {
250 J = LUSOL->iq[K];
251 I = LUSOL->ip[K];
252 L1 = LUSOL->locr[I];
253 DIAG = fabs(LUSOL->a[L1]);
254 SETMAX( DUMAX, DIAG );
255 if(DUMIN > DIAG) {
256 DUMIN = DIAG;
257 JUMIN = J;
258 }
259 }
260 }
261 else {
262 /* --------------------------------------------------------------
263 keepLU = 0.
264 Only diag(U) is stored. Set w(*) accordingly.
265 Find DUmax and DUmin, the extreme diagonals of U.
266 -------------------------------------------------------------- */
267 LDIAGU = LENA2 - LUSOL->n;
268 for(K = 1; K <= NRANK; K++) {
269 J = LUSOL->iq[K];
270 DIAG = fabs( LUSOL->a[LDIAGU + J] ); /* are in natural order */
271 LUSOL->w[J] = DIAG;
272 SETMAX( DUMAX, DIAG );
273 if(DUMIN > DIAG) {
274 DUMIN = DIAG;
275 JUMIN = J;
276 }
277 }
278 }
279 /* --------------------------------------------------------------
280 Negate w(j) if the corresponding diagonal of U is
281 too small in absolute terms or relative to the other elements
282 in the same column of U.
283
284 23 Apr 2004: TRP ensures that diags are NOT small relative to
285 other elements in their own column.
286 Much better, we can compare all diags to DUmax.
287 -------------------------------------------------------------- */
288 if((MODE == 1) && TRP) {
289 SETMAX( UTOL1, UTOL2*DUMAX );
290 }
291
292 if(KEEPLU) {
293 for(K = 1; K <= LUSOL->n; K++) {
294 J = LUSOL->iq[K];
295 if(K>NRANK)
296 DIAG = ZERO;
297 else {
298 I = LUSOL->ip[K];
299 L1 = LUSOL->locr[I];
300 DIAG = fabs(LUSOL->a[L1]);
301 }
302 if((DIAG<=UTOL1) || (DIAG<=UTOL2*LUSOL->w[J])) {
303 LUSOL_addSingularity(LUSOL, J, INFORM);
304 LUSOL->w[J] = -LUSOL->w[J];
305 }
306 }
307 }
308 else { /* keepLU = FALSE */
309 for(K = 1; K <= LUSOL->n; K++) {
310 J = LUSOL->iq[K];
311 DIAG = LUSOL->w[J];
312 if(DIAG<=UTOL1) {
313 LUSOL_addSingularity(LUSOL, J, INFORM);
314 LUSOL->w[J] = -LUSOL->w[J];
315 }
316 }
317 }
318 /* -----------------------------------------------------------------
319 Set output parameters.
320 ----------------------------------------------------------------- */
321 if(JUMIN==0)
322 DUMIN = ZERO;
323 LUSOL->luparm[LUSOL_IP_COLINDEX_DUMIN] = JUMIN;
324 LUSOL->parmlu[LUSOL_RP_MAXELEM_DIAGU] = DUMAX;
325 LUSOL->parmlu[LUSOL_RP_MINELEM_DIAGU] = DUMIN;
326 /* The matrix has been judged singular. */
327 if(LUSOL->luparm[LUSOL_IP_SINGULARITIES]>0) {
328 *INFORM = LUSOL_INFORM_LUSINGULAR;
329 NDEFIC = LUSOL->n-NRANK;
330 if((LUSOL->outstream!=NULL) && (LPRINT>=LUSOL_MSG_SINGULARITY)) {
331 LUSOL_report(LUSOL, 0, "Singular(m%cn) rank:%9d n-rank:%8d nsing:%9d\n",
332 relationChar(LUSOL->m, LUSOL->n),NRANK,NDEFIC,
333 LUSOL->luparm[LUSOL_IP_SINGULARITIES]);
334 }
335 }
336 /* Exit. */
337 LUSOL->luparm[LUSOL_IP_INFORM] = *INFORM;
338 }
339 #endif
340
341
342 /* ------------------------------------------------------------------
343 Include routines for row-based L0.
344 20 Apr 2005 Current version - KE.
345 ------------------------------------------------------------------ */
346 #include "lusol6l0.h"
347
348
349 /* ------------------------------------------------------------------
350 lu6L solves L v = v(input).
351 ------------------------------------------------------------------
352 15 Dec 2002: First version derived from lu6sol.
353 15 Dec 2002: Current version.
354 ------------------------------------------------------------------ */
LU6L(LUSOLrec * LUSOL,int * INFORM,REAL V[],int NZidx[])355 void LU6L(LUSOLrec *LUSOL, int *INFORM, REAL V[], int NZidx[])
356 {
357 int JPIV, K, L, L1, LEN, LENL, LENL0, NUML, NUML0;
358 REAL SMALL;
359 register REAL VPIV;
360 #ifdef LUSOLFastSolve
361 REAL *aptr;
362 int *iptr, *jptr;
363 #else
364 int I, J;
365 #endif
366
367 NUML0 = LUSOL->luparm[LUSOL_IP_COLCOUNT_L0];
368 LENL0 = LUSOL->luparm[LUSOL_IP_NONZEROS_L0];
369 LENL = LUSOL->luparm[LUSOL_IP_NONZEROS_L];
370 SMALL = LUSOL->parmlu[LUSOL_RP_ZEROTOLERANCE];
371 *INFORM = LUSOL_INFORM_LUSUCCESS;
372 L1 = LUSOL->lena+1;
373 for(K = 1; K <= NUML0; K++) {
374 LEN = LUSOL->lenc[K];
375 L = L1;
376 L1 -= LEN;
377 JPIV = LUSOL->indr[L1];
378 VPIV = V[JPIV];
379 if(fabs(VPIV)>SMALL) {
380 /* ***** This loop could be coded specially. */
381 #ifdef LUSOLFastSolve
382 L--;
383 for(aptr = LUSOL->a+L, iptr = LUSOL->indc+L;
384 LEN > 0; LEN--, aptr--, iptr--)
385 V[*iptr] += (*aptr) * VPIV;
386 #else
387 for(; LEN > 0; LEN--) {
388 L--;
389 I = LUSOL->indc[L];
390 V[I] += LUSOL->a[L]*VPIV;
391 }
392 #endif
393 }
394 #ifdef SetSmallToZero
395 else
396 V[JPIV] = 0;
397 #endif
398 }
399 L = (LUSOL->lena-LENL0)+1;
400 NUML = LENL-LENL0;
401 /* ***** This loop could be coded specially. */
402 #ifdef LUSOLFastSolve
403 L--;
404 for(aptr = LUSOL->a+L, jptr = LUSOL->indr+L, iptr = LUSOL->indc+L;
405 NUML > 0; NUML--, aptr--, jptr--, iptr--) {
406 if(fabs(V[*jptr])>SMALL)
407 V[*iptr] += (*aptr) * V[*jptr];
408 #ifdef SetSmallToZero
409 else
410 V[*jptr] = 0;
411 #endif
412 }
413 #else
414 for(; NUML > 0; NUML--) {
415 L--;
416 J = LUSOL->indr[L];
417 if(fabs(V[J])>SMALL) {
418 I = LUSOL->indc[L];
419 V[I] += LUSOL->a[L]*V[J];
420 }
421 #ifdef SetSmallToZero
422 else
423 V[J] = 0;
424 #endif
425 }
426 #endif
427 /* Exit. */
428 LUSOL->luparm[LUSOL_IP_INFORM] = *INFORM;
429 }
430
431 /* ==================================================================
432 lu6LD assumes lu1fac has computed factors A = LU of a
433 symmetric definite or quasi-definite matrix A,
434 using Threshold Symmetric Pivoting (TSP), luparm(6) = 3,
435 or Threshold Diagonal Pivoting (TDP), luparm(6) = 4.
436 It also assumes that no updates have been performed.
437 In such cases, U = D L', where D = diag(U).
438 lu6LDL returns v as follows:
439
440 mode
441 1 v solves L D v = v(input).
442 2 v solves L|D|v = v(input).
443 ------------------------------------------------------------------
444 15 Dec 2002: First version of lu6LD.
445 15 Dec 2002: Current version.
446 ================================================================== */
LU6LD(LUSOLrec * LUSOL,int * INFORM,int MODE,REAL V[],int NZidx[])447 void LU6LD(LUSOLrec *LUSOL, int *INFORM, int MODE, REAL V[], int NZidx[])
448 {
449 int IPIV, K, L, L1, LEN, NUML0;
450 REAL DIAG, SMALL;
451 register REAL VPIV;
452 #ifdef LUSOLFastSolve
453 REAL *aptr;
454 int *jptr;
455 #else
456 int J;
457 #endif
458
459 /* Solve L D v(new) = v or L|D|v(new) = v, depending on mode.
460 The code for L is the same as in lu6L,
461 but when a nonzero entry of v arises, we divide by
462 the corresponding entry of D or |D|. */
463 NUML0 = LUSOL->luparm[LUSOL_IP_COLCOUNT_L0];
464 SMALL = LUSOL->parmlu[LUSOL_RP_ZEROTOLERANCE];
465 *INFORM = LUSOL_INFORM_LUSUCCESS;
466 L1 = LUSOL->lena+1;
467 for(K = 1; K <= NUML0; K++) {
468 LEN = LUSOL->lenc[K];
469 L = L1;
470 L1 -= LEN;
471 IPIV = LUSOL->indr[L1];
472 VPIV = V[IPIV];
473 if(fabs(VPIV)>SMALL) {
474 /* ***** This loop could be coded specially. */
475 #ifdef LUSOLFastSolve
476 L--;
477 for(aptr = LUSOL->a+L, jptr = LUSOL->indc+L;
478 LEN > 0; LEN--, aptr--, jptr--)
479 V[*jptr] += (*aptr)*VPIV;
480 #else
481 for(; LEN > 0; LEN--) {
482 L--;
483 J = LUSOL->indc[L];
484 V[J] += LUSOL->a[L]*VPIV;
485 }
486 #endif
487 /* Find diag = U(ipiv,ipiv) and divide by diag or |diag|. */
488 L = LUSOL->locr[IPIV];
489 DIAG = LUSOL->a[L];
490 if(MODE==2)
491 DIAG = fabs(DIAG);
492 V[IPIV] = VPIV/DIAG;
493 }
494 #ifdef SetSmallToZero
495 else
496 V[IPIV] = 0;
497 #endif
498 }
499 }
500
501
502 /* ==================================================================
503 lu6Lt solves L'v = v(input).
504 ------------------------------------------------------------------
505 15 Dec 2002: First version derived from lu6sol.
506 15 Dec 2002: Current version.
507 ================================================================== */
LU6LT(LUSOLrec * LUSOL,int * INFORM,REAL V[],int NZidx[])508 void LU6LT(LUSOLrec *LUSOL, int *INFORM, REAL V[], int NZidx[])
509 {
510 #ifdef DoTraceL0
511 REAL TEMP;
512 #endif
513 int K, L, L1, L2, LEN, LENL, LENL0, NUML0;
514 REAL SMALL;
515 register REALXP SUM;
516 register REAL HOLD;
517 #if (defined LUSOLFastSolve) && !(defined DoTraceL0)
518 REAL *aptr;
519 int *iptr, *jptr;
520 #else
521 int I, J;
522 #endif
523
524 NUML0 = LUSOL->luparm[LUSOL_IP_COLCOUNT_L0];
525 LENL0 = LUSOL->luparm[LUSOL_IP_NONZEROS_L0];
526 LENL = LUSOL->luparm[LUSOL_IP_NONZEROS_L];
527 SMALL = LUSOL->parmlu[LUSOL_RP_ZEROTOLERANCE];
528 *INFORM = LUSOL_INFORM_LUSUCCESS;
529 L1 = (LUSOL->lena-LENL)+1;
530 L2 = LUSOL->lena-LENL0;
531
532 /* ***** This loop could be coded specially. */
533 #if (defined LUSOLFastSolve) && !(defined DoTraceL0)
534 for(L = L1, aptr = LUSOL->a+L1, iptr = LUSOL->indr+L1, jptr = LUSOL->indc+L1;
535 L <= L2; L++, aptr++, iptr++, jptr++) {
536 HOLD = V[*jptr];
537 if(fabs(HOLD)>SMALL)
538 V[*iptr] += (*aptr)*HOLD;
539 #ifdef SetSmallToZero
540 else
541 V[*jptr] = 0;
542 #endif
543 }
544 #else
545 for(L = L1; L <= L2; L++) {
546 J = LUSOL->indc[L];
547 HOLD = V[J];
548 if(fabs(HOLD)>SMALL) {
549 I = LUSOL->indr[L];
550 V[I] += LUSOL->a[L]*HOLD;
551 }
552 #ifdef SetSmallToZero
553 else
554 V[J] = 0;
555 #endif
556 }
557 #endif
558
559 /* Do row-based L0 version, if available */
560 if((LUSOL->L0 != NULL) ||
561 ((LUSOL->luparm[LUSOL_IP_BTRANCOUNT] == 0) && LU1L0(LUSOL, &(LUSOL->L0), INFORM))) {
562 LU6L0T_v(LUSOL, LUSOL->L0, V, NZidx, INFORM);
563 }
564
565 /* Alternatively, do the standard column-based L0 version */
566 else {
567 /* Perform loop over columns */
568 for(K = NUML0; K >= 1; K--) {
569 SUM = ZERO;
570 LEN = LUSOL->lenc[K];
571 L1 = L2+1;
572 L2 += LEN;
573 /* ***** This loop could be coded specially. */
574 #if (defined LUSOLFastSolve) && !(defined DoTraceL0)
575 for(L = L1, aptr = LUSOL->a+L1, jptr = LUSOL->indc+L1;
576 L <= L2; L++, aptr++, jptr++)
577 SUM += (*aptr) * V[*jptr];
578 #else
579 for(L = L1; L <= L2; L++) {
580 J = LUSOL->indc[L];
581 #ifndef DoTraceL0
582 SUM += LUSOL->a[L]*V[J];
583 #else
584 TEMP = V[LUSOL->indr[L1]] + SUM;
585 SUM += LUSOL->a[L]*V[J];
586 printf("V[%3d] = V[%3d] + L[%d,%d]*V[%3d]\n", LUSOL->indr[L1], LUSOL->indr[L1], J,LUSOL->indr[L1], J);
587 printf("%6g = %6g + %6g*%6g\n", V[LUSOL->indr[L1]] + SUM, TEMP, LUSOL->a[L], V[J]);
588 #endif
589 }
590 #endif
591 V[LUSOL->indr[L1]] += (REAL) SUM;
592 }
593 }
594
595 /* Exit. */
596 LUSOL->luparm[LUSOL_IP_INFORM] = *INFORM;
597 }
598
print_L0(LUSOLrec * LUSOL)599 void print_L0(LUSOLrec *LUSOL)
600 {
601 int I, J, K, L, L1, L2, LEN, LENL0, NUML0;
602 REAL *denseL0 = (REAL*) calloc(LUSOL->m+1, (LUSOL->n+1)*sizeof(*denseL0));
603
604 NUML0 = LUSOL->luparm[LUSOL_IP_COLCOUNT_L0];
605 LENL0 = LUSOL->luparm[LUSOL_IP_NONZEROS_L0];
606
607 L2 = LUSOL->lena-LENL0;
608 for(K = NUML0; K >= 1; K--) {
609 LEN = LUSOL->lenc[K];
610 L1 = L2+1;
611 L2 += LEN;
612 for(L = L1; L <= L2; L++) {
613 I = LUSOL->indc[L];
614 I = LUSOL->ipinv[I]; /* Undo row mapping */
615 J = LUSOL->indr[L];
616 denseL0[(LUSOL->n+1)*(J-1) + I] = LUSOL->a[L];
617 }
618 }
619
620 /* for(I = 1; I <= LUSOL->n; I++) {
621 ** for(J = 1; J <= LUSOL->m; J++)
622 ** fprintf(stdout, "%10g", denseL0[(LUSOL->n+1)*(J-1) + I]);
623 ** fprintf(stdout, "\n");
624 **} */
625 LUSOL_FREE(denseL0);
626 }
627
628
629 /* ------------------------------------------------------------------
630 Include routines for column-based U.
631 5 Feb 2006 Current version - KE.
632 ------------------------------------------------------------------ */
633 #include "lusol6u.h"
634
635
636 /* ==================================================================
637 lu6U solves U w = v. v is not altered.
638 ------------------------------------------------------------------
639 15 Dec 2002: First version derived from lu6sol.
640 15 Dec 2002: Current version.
641 ================================================================== */
LU6U(LUSOLrec * LUSOL,int * INFORM,REAL V[],REAL W[],int NZidx[])642 void LU6U(LUSOLrec *LUSOL, int *INFORM, REAL V[], REAL W[], int NZidx[])
643 {
644 /* Do column-based U version, if available */
645 if((LUSOL->U != NULL) ||
646 ((LUSOL->luparm[LUSOL_IP_FTRANCOUNT] == 0) && LU1U0(LUSOL, &(LUSOL->U), INFORM))) {
647 LU6U0_v(LUSOL, LUSOL->U, V, W, NZidx, INFORM);
648 }
649 /* Alternatively, do the standard column-based L0 version */
650 else {
651 int I, J, K, KLAST, L, L1, L2, L3, NRANK, NRANK1;
652 REAL SMALL;
653 register REALXP T;
654 #ifdef LUSOLFastSolve
655 REAL *aptr;
656 int *jptr;
657 #endif
658 NRANK = LUSOL->luparm[LUSOL_IP_RANK_U];
659 SMALL = LUSOL->parmlu[LUSOL_RP_ZEROTOLERANCE];
660 *INFORM = LUSOL_INFORM_LUSUCCESS;
661 NRANK1 = NRANK+1;
662 /* Find the first nonzero in v(1:nrank), counting backwards. */
663 for(KLAST = NRANK; KLAST >= 1; KLAST--) {
664 I = LUSOL->ip[KLAST];
665 if(fabs(V[I])>SMALL)
666 break;
667 }
668 L = LUSOL->n;
669 #ifdef LUSOLFastSolve
670 for(K = KLAST+1, jptr = LUSOL->iq+K; K <= L; K++, jptr++)
671 W[*jptr] = ZERO;
672 #else
673 for(K = KLAST+1; K <= L; K++) {
674 J = LUSOL->iq[K];
675 W[J] = ZERO;
676 }
677 #endif
678 /* Do the back-substitution, using rows 1:klast of U. */
679 for(K = KLAST; K >= 1; K--) {
680 I = LUSOL->ip[K];
681 T = V[I];
682 L1 = LUSOL->locr[I];
683 L2 = L1+1;
684 L3 = (L1+LUSOL->lenr[I])-1;
685 /* ***** This loop could be coded specially. */
686 #ifdef LUSOLFastSolve
687 for(L = L2, aptr = LUSOL->a+L2, jptr = LUSOL->indr+L2;
688 L <= L3; L++, aptr++, jptr++)
689 T -= (*aptr) * W[*jptr];
690 #else
691 for(L = L2; L <= L3; L++) {
692 J = LUSOL->indr[L];
693 T -= LUSOL->a[L]*W[J];
694 }
695 #endif
696 J = LUSOL->iq[K];
697 if(fabs((REAL) T)<=SMALL)
698 T = ZERO;
699 else
700 T /= LUSOL->a[L1];
701 W[J] = (REAL) T;
702 }
703 /* Compute residual for overdetermined systems. */
704 T = ZERO;
705 for(K = NRANK1; K <= LUSOL->m; K++) {
706 I = LUSOL->ip[K];
707 T += fabs(V[I]);
708 }
709 /* Exit. */
710 if(T>ZERO)
711 *INFORM = LUSOL_INFORM_LUSINGULAR;
712 LUSOL->luparm[LUSOL_IP_INFORM] = *INFORM;
713 LUSOL->parmlu[LUSOL_RP_RESIDUAL_U] = (REAL) T;
714 }
715 }
716
717 /* ==================================================================
718 lu6Ut solves U'v = w. w is destroyed.
719 ------------------------------------------------------------------
720 15 Dec 2002: First version derived from lu6sol.
721 15 Dec 2002: Current version.
722 ================================================================== */
LU6UT(LUSOLrec * LUSOL,int * INFORM,REAL V[],REAL W[],int NZidx[])723 void LU6UT(LUSOLrec *LUSOL, int *INFORM, REAL V[], REAL W[], int NZidx[])
724 {
725 int I, J, K, L, L1, L2, NRANK, NRANK1,
726 *ip = LUSOL->ip + 1, *iq = LUSOL->iq + 1;
727 REAL SMALL;
728 register REAL T;
729 #ifdef LUSOLFastSolve
730 REAL *aptr;
731 int *jptr;
732 #endif
733
734 NRANK = LUSOL->luparm[LUSOL_IP_RANK_U];
735 SMALL = LUSOL->parmlu[LUSOL_RP_ZEROTOLERANCE];
736 *INFORM = LUSOL_INFORM_LUSUCCESS;
737 NRANK1 = NRANK+1;
738 L = LUSOL->m;
739 #ifdef LUSOLFastSolve
740 for(K = NRANK1, jptr = LUSOL->ip+K; K <= L; K++, jptr++)
741 V[*jptr] = ZERO;
742 #else
743 for(K = NRANK1; K <= L; K++) {
744 I = LUSOL->ip[K];
745 V[I] = ZERO;
746 }
747 #endif
748 /* Do the forward-substitution, skipping columns of U(transpose)
749 when the associated element of w(*) is negligible. */
750 #if 0
751 for(K = 1; K <= NRANK; K++) {
752 I = LUSOL->ip[K];
753 J = LUSOL->iq[K];
754 #else
755 for(K = 1; K <= NRANK; K++, ip++, iq++) {
756 I = *ip;
757 J = *iq;
758 #endif
759 T = W[J];
760 if(fabs(T)<=SMALL) {
761 V[I] = ZERO;
762 continue;
763 }
764 L1 = LUSOL->locr[I];
765 T /= LUSOL->a[L1];
766 V[I] = T;
767 L2 = (L1+LUSOL->lenr[I])-1;
768 L1++;
769 /* ***** This loop could be coded specially. */
770 #ifdef LUSOLFastSolve
771 for(L = L1, aptr = LUSOL->a+L1, jptr = LUSOL->indr+L1;
772 L <= L2; L++, aptr++, jptr++)
773 W[*jptr] -= T * (*aptr);
774 #else
775 for(L = L1; L <= L2; L++) {
776 J = LUSOL->indr[L];
777 W[J] -= T*LUSOL->a[L];
778 }
779 #endif
780 }
781 /* Compute residual for overdetermined systems. */
782 T = ZERO;
783 for(K = NRANK1; K <= LUSOL->n; K++) {
784 J = LUSOL->iq[K];
785 T += fabs(W[J]);
786 }
787 /* Exit. */
788 if(T>ZERO)
789 *INFORM = LUSOL_INFORM_LUSINGULAR;
790 LUSOL->luparm[LUSOL_IP_INFORM] = *INFORM;
791 LUSOL->parmlu[LUSOL_RP_RESIDUAL_U] = T;
792 }
793
794 /* ==================================================================
795 lu6sol uses the factorization A = L U as follows:
796 ------------------------------------------------------------------
797 mode
798 1 v solves L v = v(input). w is not touched.
799 2 v solves L'v = v(input). w is not touched.
800 3 w solves U w = v. v is not altered.
801 4 v solves U'v = w. w is destroyed.
802 5 w solves A w = v. v is altered as in 1.
803 6 v solves A'v = w. w is destroyed.
804
805 If mode = 3,4,5,6, v and w must not be the same arrays.
806 If lu1fac has just been used to factorize a symmetric matrix A
807 (which must be definite or quasi-definite), the factors A = L U
808 may be regarded as A = LDL', where D = diag(U). In such cases,
809
810 mode
811 7 v solves A v = L D L'v = v(input). w is not touched.
812 8 v solves L |D| L'v = v(input). w is not touched.
813
814 ip(*), iq(*) hold row and column numbers in pivotal order.
815 lenc(k) is the length of the k-th column of initial L.
816 lenr(i) is the length of the i-th row of U.
817 locc(*) is not used.
818 locr(i) is the start of the i-th row of U.
819
820 U is assumed to be in upper-trapezoidal form (nrank by n).
821 The first entry for each row is the diagonal element
822 (according to the permutations ip, iq). It is stored at
823 location locr(i) in a(*), indr(*).
824
825 On exit, inform = 0 except as follows.
826 if(mode = 3,4,5,6 and if U (and hence A) is singular,)
827 inform = 1 if there is a nonzero residual in solving the system
828 involving U. parmlu(20) returns the norm of the residual.
829 ------------------------------------------------------------------
830 July 1987: Early version.
831 09 May 1988: f77 version.
832 27 Apr 2000: Abolished the dreaded "computed go to".
833 But hard to change other "go to"s to "if then else".
834 15 Dec 2002: lu6L, lu6Lt, lu6U, lu6Ut added to modularize lu6sol.
835 ================================================================== */
836 void LU6SOL(LUSOLrec *LUSOL, int MODE, REAL V[], REAL W[], int NZidx[], int *INFORM)
837 {
838 if(MODE==LUSOL_SOLVE_Lv_v) { /* Solve L v(new) = v. */
839 LU6L(LUSOL, INFORM,V, NZidx);
840 }
841 else if(MODE==LUSOL_SOLVE_Ltv_v) { /* Solve L'v(new) = v. */
842 LU6LT(LUSOL, INFORM,V, NZidx);
843 }
844 else if(MODE==LUSOL_SOLVE_Uw_v) { /* Solve U w = v. */
845 LU6U(LUSOL, INFORM,V,W, NZidx);
846 }
847 else if(MODE==LUSOL_SOLVE_Utv_w) { /* Solve U'v = w. */
848 LU6UT(LUSOL, INFORM,V,W, NZidx);
849 }
850 else if(MODE==LUSOL_SOLVE_Aw_v) { /* Solve A w = v (i.e. FTRAN) */
851 LU6L(LUSOL, INFORM,V, NZidx); /* via L v(new) = v */
852 LU6U(LUSOL, INFORM,V,W, NULL); /* ... and U w = v(new). */
853 }
854 else if(MODE==LUSOL_SOLVE_Atv_w) { /* Solve A'v = w (i.e. BTRAN) */
855 LU6UT(LUSOL, INFORM,V,W, NZidx); /* via U'v = w */
856 LU6LT(LUSOL, INFORM,V, NULL); /* ... and L'v(new) = v. */
857 }
858 else if(MODE==LUSOL_SOLVE_Av_v) { /* Solve LDv(bar) = v */
859 LU6LD(LUSOL, INFORM,1,V, NZidx); /* and L'v(new) = v(bar). */
860 LU6LT(LUSOL, INFORM,V, NULL);
861 }
862 else if(MODE==LUSOL_SOLVE_LDLtv_v) { /* Solve L|D|v(bar) = v */
863 LU6LD(LUSOL, INFORM,2,V, NZidx); /* and L'v(new) = v(bar). */
864 LU6LT(LUSOL, INFORM,V, NULL);
865 }
866 }
867
868