1
2 /* ==================================================================
3 lu1DCP factors a dense m x n matrix A by Gaussian elimination,
4 using Complete Pivoting (row and column interchanges) for stability.
5 This version also uses column interchanges if all elements in a
6 pivot column are smaller than (or equal to) "small". Such columns
7 are changed to zero and permuted to the right-hand end.
8 As in LINPACK's dgefa, ipvt(!) keeps track of pivot rows.
9 Rows of U are interchanged, but we don't have to physically
10 permute rows of L. In contrast, column interchanges are applied
11 directly to the columns of both L and U, and to the column
12 permutation vector iq(*).
13 ------------------------------------------------------------------
14 On entry:
15 a Array holding the matrix A to be factored.
16 lda The leading dimension of the array a.
17 m The number of rows in A.
18 n The number of columns in A.
19 small A drop tolerance. Must be zero or positive.
20
21 On exit:
22 a An upper triangular matrix and the multipliers
23 which were used to obtain it.
24 The factorization can be written A = L*U where
25 L is a product of permutation and unit lower
26 triangular matrices and U is upper triangular.
27 nsing Number of singularities detected.
28 ipvt Records the pivot rows.
29 iq A vector to which column interchanges are applied.
30 ------------------------------------------------------------------
31 01 May 2002: First dense Complete Pivoting, derived from lu1DPP.
32 07 May 2002: Another break needed at end of first loop.
33 07 May 2002: Current version of lu1DCP.
34 ================================================================== */
LU1DCP(LUSOLrec * LUSOL,REAL DA[],int LDA,int M,int N,REAL SMALL,int * NSING,int IPVT[],int IX[])35 void LU1DCP(LUSOLrec *LUSOL, REAL DA[], int LDA, int M, int N, REAL SMALL,
36 int *NSING, int IPVT[], int IX[])
37 {
38
39 int I, J, K, KP1, L, LAST, LENCOL, IMAX, JMAX, JLAST, JNEW;
40 REAL AIJMAX, AJMAX;
41 register REAL T;
42 #ifdef LUSOLFastDenseIndex
43 register REAL *DA1, *DA2;
44 int IDA1, IDA2;
45 #else
46 register int IDA1, IDA2;
47 #endif
48
49 *NSING = 0;
50 LENCOL = M+1;
51 LAST = N;
52 /* -----------------------------------------------------------------
53 Start of elimination loop.
54 ----------------------------------------------------------------- */
55 for(K = 1; K <= N; K++) {
56 KP1 = K+1;
57 LENCOL--;
58 /* Find the biggest aij in row imax and column jmax. */
59 AIJMAX = ZERO;
60 IMAX = K;
61 JMAX = K;
62 JLAST = LAST;
63 for(J = K; J <= JLAST; J++) {
64 x10:
65 L = idamaxlpsolve(LENCOL,DA+DAPOS(K,J)-LUSOL_ARRAYOFFSET,1)+K-1;
66 AJMAX = fabs(DA[DAPOS(L,J)]);
67 if(AJMAX<=SMALL) {
68 /* ========================================================
69 Do column interchange, changing old column to zero.
70 Reduce "last" and try again with same j.
71 ======================================================== */
72 (*NSING)++;
73 JNEW = IX[LAST];
74 IX[LAST] = IX[J];
75 IX[J] = JNEW;
76 #ifdef LUSOLFastDenseIndex
77 DA1 = DA+DAPOS(0,LAST);
78 DA2 = DA+DAPOS(0,J);
79 for(I = 1; I <= K-1; I++) {
80 DA1++;
81 DA2++;
82 T = *DA1;
83 *DA1 = *DA2;
84 *DA2 = T;
85 #else
86 for(I = 1; I <= K-1; I++) {
87 IDA1 = DAPOS(I,LAST);
88 IDA2 = DAPOS(I,J);
89 T = DA[IDA1];
90 DA[IDA1] = DA[IDA2];
91 DA[IDA2] = T;
92 #endif
93 }
94 #ifdef LUSOLFastDenseIndex
95 for(I = K; I <= M; I++) {
96 DA1++;
97 DA2++;
98 T = *DA1;
99 *DA1 = ZERO;
100 *DA2 = T;
101 #else
102 for(I = K; I <= M; I++) {
103 IDA1 = DAPOS(I,LAST);
104 IDA2 = DAPOS(I,J);
105 T = DA[IDA1];
106 DA[IDA1] = ZERO;
107 DA[IDA2] = T;
108 #endif
109 }
110 LAST--;
111 if(J<=LAST)
112 goto x10;
113 break;
114 }
115 /* Check if this column has biggest aij so far. */
116 if(AIJMAX<AJMAX) {
117 AIJMAX = AJMAX;
118 IMAX = L;
119 JMAX = J;
120 }
121 if(J>=LAST)
122 break;
123 }
124 IPVT[K] = IMAX;
125 if(JMAX!=K) {
126 /* ==========================================================
127 Do column interchange (k and jmax).
128 ========================================================== */
129 JNEW = IX[JMAX];
130 IX[JMAX] = IX[K];
131 IX[K] = JNEW;
132 #ifdef LUSOLFastDenseIndex
133 DA1 = DA+DAPOS(0,JMAX);
134 DA2 = DA+DAPOS(0,K);
135 for(I = 1; I <= M; I++) {
136 DA1++;
137 DA2++;
138 T = *DA1;
139 *DA1 = *DA2;
140 *DA2 = T;
141 #else
142 for(I = 1; I <= M; I++) {
143 IDA1 = DAPOS(I,JMAX);
144 IDA2 = DAPOS(I,K);
145 T = DA[IDA1];
146 DA[IDA1] = DA[IDA2];
147 DA[IDA2] = T;
148 #endif
149 }
150 }
151 if(M>K) {
152 /* ===========================================================
153 Do row interchange if necessary.
154 =========================================================== */
155 if(IMAX!=K) {
156 IDA1 = DAPOS(IMAX,K);
157 IDA2 = DAPOS(K,K);
158 T = DA[IDA1];
159 DA[IDA1] = DA[IDA2];
160 DA[IDA2] = T;
161 }
162 /* ===========================================================
163 Compute multipliers.
164 Do row elimination with column indexing.
165 =========================================================== */
166 T = -ONE/DA[DAPOS(K,K)];
167 dscallpsolve(M-K,T,DA+DAPOS(KP1,K)-LUSOL_ARRAYOFFSET,1);
168 for(J = KP1; J <= LAST; J++) {
169 IDA1 = DAPOS(IMAX,J);
170 T = DA[IDA1];
171 if(IMAX!=K) {
172 IDA2 = DAPOS(K,J);
173 DA[IDA1] = DA[IDA2];
174 DA[IDA2] = T;
175 }
176 daxpylpsolve(M-K,T,DA+DAPOS(KP1,K)-LUSOL_ARRAYOFFSET,1,
177 DA+DAPOS(KP1,J)-LUSOL_ARRAYOFFSET,1);
178 }
179 }
180 else
181 break;
182 if(K>=LAST)
183 break;
184 }
185 /* Set ipvt(*) for singular rows. */
186 for(K = LAST+1; K <= M; K++)
187 IPVT[K] = K;
188
189 }
190
191 /* ==================================================================
192 lu1DPP factors a dense m x n matrix A by Gaussian elimination,
193 using row interchanges for stability, as in dgefa from LINPACK.
194 This version also uses column interchanges if all elements in a
195 pivot column are smaller than (or equal to) "small". Such columns
196 are changed to zero and permuted to the right-hand end.
197 As in LINPACK, ipvt(*) keeps track of pivot rows.
198 Rows of U are interchanged, but we don't have to physically
199 permute rows of L. In contrast, column interchanges are applied
200 directly to the columns of both L and U, and to the column
201 permutation vector iq(*).
202 ------------------------------------------------------------------
203 On entry:
204 a Array holding the matrix A to be factored.
205 lda The leading dimension of the array a.
206 m The number of rows in A.
207 n The number of columns in A.
208 small A drop tolerance. Must be zero or positive.
209
210 On exit:
211 a An upper triangular matrix and the multipliers
212 which were used to obtain it.
213 The factorization can be written A = L*U where
214 L is a product of permutation and unit lower
215 triangular matrices and U is upper triangular.
216 nsing Number of singularities detected.
217 ipvt Records the pivot rows.
218 iq A vector to which column interchanges are applied.
219 ------------------------------------------------------------------
220 02 May 1989: First version derived from dgefa
221 in LINPACK (version dated 08/14/78).
222 05 Feb 1994: Generalized to treat rectangular matrices
223 and use column interchanges when necessary.
224 ipvt is retained, but column permutations are applied
225 directly to iq(*).
226 21 Dec 1994: Bug found via example from Steve Dirkse.
227 Loop 100 added to set ipvt(*) for singular rows.
228 ================================================================== */
229 void LU1DPP(LUSOLrec *LUSOL, REAL DA[], int LDA, int M, int N, REAL SMALL,
230 int *NSING, int IPVT[], int IX[])
231 {
232 int I, J, K, KP1, L, LAST, LENCOL;
233 register REAL T;
234 #ifdef LUSOLFastDenseIndex
235 register REAL *DA1, *DA2;
236 int IDA1, IDA2;
237 #else
238 register int IDA1, IDA2;
239 #endif
240
241 *NSING = 0;
242 K = 1;
243 LAST = N;
244 /* ------------------------------------------------------------------
245 Start of elimination loop.
246 ------------------------------------------------------------------ */
247 x10:
248 KP1 = K+1;
249 LENCOL = (M-K)+1;
250 /* Find l, the pivot row. */
251 L = (idamaxlpsolve(LENCOL,DA+DAPOS(K,K)-LUSOL_ARRAYOFFSET,1)+K)-1;
252 IPVT[K] = L;
253 if(fabs(DA[DAPOS(L,K)])<=SMALL) {
254 /* ===============================================================
255 Do column interchange, changing old pivot column to zero.
256 Reduce "last" and try again with same k.
257 =============================================================== */
258 (*NSING)++;
259 J = IX[LAST];
260 IX[LAST] = IX[K];
261 IX[K] = J;
262 #ifdef LUSOLFastDenseIndex
263 DA1 = DA+DAPOS(0,LAST);
264 DA2 = DA+DAPOS(0,K);
265 for(I = 1; I <= K-1; I++) {
266 DA1++;
267 DA2++;
268 T = *DA1;
269 *DA1 = *DA2;
270 *DA2 = T;
271 #else
272 for(I = 1; I <= K-1; I++) {
273 IDA1 = DAPOS(I,LAST);
274 IDA2 = DAPOS(I,K);
275 T = DA[IDA1];
276 DA[IDA1] = DA[IDA2];
277 DA[IDA2] = T;
278 #endif
279 }
280 #ifdef LUSOLFastDenseIndex
281 for(I = K; I <= M; I++) {
282 DA1++;
283 DA2++;
284 T = *DA1;
285 *DA1 = ZERO;
286 *DA2 = T;
287 #else
288 for(I = K; I <= M; I++) {
289 IDA1 = DAPOS(I,LAST);
290 IDA2 = DAPOS(I,K);
291 T = DA[IDA1];
292 DA[IDA1] = ZERO;
293 DA[IDA2] = T;
294 #endif
295 }
296 LAST = LAST-1;
297 if(K<=LAST)
298 goto x10;
299 }
300 else if(M>K) {
301 /* ===============================================================
302 Do row interchange if necessary.
303 =============================================================== */
304 if(L!=K) {
305 IDA1 = DAPOS(L,K);
306 IDA2 = DAPOS(K,K);
307 T = DA[IDA1];
308 DA[IDA1] = DA[IDA2];
309 DA[IDA2] = T;
310 }
311 /* ===============================================================
312 Compute multipliers.
313 Do row elimination with column indexing.
314 =============================================================== */
315 T = -ONE/DA[DAPOS(K,K)];
316 dscallpsolve(M-K,T,DA+DAPOS(KP1,K)-LUSOL_ARRAYOFFSET,1);
317 for(J = KP1; J <= LAST; J++) {
318 IDA1 = DAPOS(L,J);
319 T = DA[IDA1];
320 if(L!=K) {
321 IDA2 = DAPOS(K,J);
322 DA[IDA1] = DA[IDA2];
323 DA[IDA2] = T;
324 }
325 daxpylpsolve(M-K,T,DA+DAPOS(KP1,K)-LUSOL_ARRAYOFFSET,1,
326 DA+DAPOS(KP1,J)-LUSOL_ARRAYOFFSET,1);
327 }
328 K++;
329 if(K<=LAST)
330 goto x10;
331 }
332 /* Set ipvt(*) for singular rows. */
333 for(K = LAST+1; K <= M; K++)
334 IPVT[K] = K;
335
336 }
337
338
339 /* ==================================================================
340 lu1pq1 constructs a permutation iperm from the array len.
341 ------------------------------------------------------------------
342 On entry:
343 len(i) holds the number of nonzeros in the i-th row (say)
344 of an m by n matrix.
345 num(*) can be anything (workspace).
346
347 On exit:
348 iperm contains a list of row numbers in the order
349 rows of length 0, rows of length 1,..., rows of length n.
350 loc(nz) points to the first row containing nz nonzeros,
351 nz = 1, n.
352 inv(i) points to the position of row i within iperm(*).
353 ================================================================== */
354 void LU1PQ1(LUSOLrec *LUSOL, int M, int N, int LEN[],
355 int IPERM[], int LOC[], int INV[], int NUM[])
356 {
357 int NZEROS, NZ, I, L;
358
359 /* Count the number of rows of each length. */
360 NZEROS = 0;
361 for(NZ = 1; NZ <= N; NZ++) {
362 NUM[NZ] = 0;
363 LOC[NZ] = 0;
364 }
365 for(I = 1; I <= M; I++) {
366 NZ = LEN[I];
367 if(NZ==0)
368 NZEROS++;
369 else
370 NUM[NZ]++;
371 }
372 /* Set starting locations for each length. */
373 L = NZEROS+1;
374 for(NZ = 1; NZ <= N; NZ++) {
375 LOC[NZ] = L;
376 L += NUM[NZ];
377 NUM[NZ] = 0;
378 }
379 /* Form the list. */
380 NZEROS = 0;
381 for(I = 1; I <= M; I++) {
382 NZ = LEN[I];
383 if(NZ==0) {
384 NZEROS++;
385 IPERM[NZEROS] = I;
386 }
387 else {
388 L = LOC[NZ]+NUM[NZ];
389 IPERM[L] = I;
390 NUM[NZ]++;
391 }
392 }
393 /* Define the inverse of iperm. */
394 for(L = 1; L <= M; L++) {
395 I = IPERM[L];
396 INV[I] = L;
397 }
398 }
399
400 /* ==================================================================
401 lu1pq2 frees the space occupied by the pivot row,
402 and updates the column permutation iq.
403 Also used to free the pivot column and update the row perm ip.
404 ------------------------------------------------------------------
405 nzpiv (input) is the length of the pivot row (or column).
406 nzchng (output) is the net change in total nonzeros.
407 ------------------------------------------------------------------
408 14 Apr 1989 First version.
409 ================================================================== */
410 void LU1PQ2(LUSOLrec *LUSOL, int NZPIV, int *NZCHNG,
411 int IND[], int LENOLD[], int LENNEW[], int IXLOC[], int IX[], int IXINV[])
412 {
413 int LR, J, NZ, NZNEW, L, NEXT, LNEW, JNEW;
414
415 *NZCHNG = 0;
416 for(LR = 1; LR <= NZPIV; LR++) {
417 J = IND[LR];
418 IND[LR] = 0;
419 NZ = LENOLD[LR];
420 NZNEW = LENNEW[J];
421 if(NZ!=NZNEW) {
422 L = IXINV[J];
423 *NZCHNG = (*NZCHNG+NZNEW)-NZ;
424 /* l above is the position of column j in iq (so j = iq(l)). */
425 if(NZ<NZNEW) {
426 /* Column j has to move towards the end of iq. */
427 x110:
428 NEXT = NZ+1;
429 LNEW = IXLOC[NEXT]-1;
430 if(LNEW!=L) {
431 JNEW = IX[LNEW];
432 IX[L] = JNEW;
433 IXINV[JNEW] = L;
434 }
435 L = LNEW;
436 IXLOC[NEXT] = LNEW;
437 NZ = NEXT;
438 if(NZ<NZNEW)
439 goto x110;
440 }
441 else {
442 /* Column j has to move towards the front of iq. */
443 x120:
444 LNEW = IXLOC[NZ];
445 if(LNEW!=L) {
446 JNEW = IX[LNEW];
447 IX[L] = JNEW;
448 IXINV[JNEW] = L;
449 }
450 L = LNEW;
451 IXLOC[NZ] = LNEW+1;
452 NZ = NZ-1;
453 if(NZ>NZNEW)
454 goto x120;
455 }
456 IX[LNEW] = J;
457 IXINV[J] = LNEW;
458 }
459 }
460 }
461
462 /* ==================================================================
463 lu1pq3 looks at the permutation iperm(*) and moves any entries
464 to the end whose corresponding length len(*) is zero.
465 ------------------------------------------------------------------
466 09 Feb 1994: Added work array iw(*) to improve efficiency.
467 ================================================================== */
468 void LU1PQ3(LUSOLrec *LUSOL, int MN, int LEN[], int IPERM[], int IW[], int *NRANK)
469 {
470 int NZEROS, K, I;
471
472 *NRANK = 0;
473 NZEROS = 0;
474 for(K = 1; K <= MN; K++) {
475 I = IPERM[K];
476 if(LEN[I]==0) {
477 NZEROS++;
478 IW[NZEROS] = I;
479 }
480 else {
481 (*NRANK)++;
482 IPERM[*NRANK] = I;
483 }
484 }
485 for(K = 1; K <= NZEROS; K++)
486 IPERM[(*NRANK)+K] = IW[K];
487 }
488
489 /* ==================================================================
490 lu1rec
491 ------------------------------------------------------------------
492 On exit:
493 ltop is the length of useful entries in ind(*), a(*).
494 ind(ltop+1) is "i" such that len(i), loc(i) belong to the last
495 item in ind(*), a(*).
496 ------------------------------------------------------------------
497 00 Jun 1983: Original version of lu1rec followed John Reid's
498 compression routine in LA05. It recovered
499 space in ind(*) and optionally a(*)
500 by eliminating entries with ind(l) = 0.
501 The elements of ind(*) could not be negative.
502 If len(i) was positive, entry i contained
503 that many elements, starting at loc(i).
504 Otherwise, entry i was eliminated.
505 23 Mar 2001: Realised we could have len(i) = 0 in rare cases!
506 (Mostly during TCP when the pivot row contains
507 a column of length 1 that couldn't be a pivot.)
508 Revised storage scheme to
509 keep entries with ind(l) > 0,
510 squeeze out entries with -n <= ind(l) <= 0,
511 and to allow len(i) = 0.
512 Empty items are moved to the end of the compressed
513 ind(*) and/or a(*) arrays are given one empty space.
514 Items with len(i) < 0 are still eliminated.
515 27 Mar 2001: Decided to use only ind(l) > 0 and = 0 in lu1fad.
516 Still have to keep entries with len(i) = 0.
517 ================================================================== */
518 void LU1REC(LUSOLrec *LUSOL, int N, MYBOOL REALS, int *LTOP,
519 int IND[], int LEN[], int LOC[])
520 {
521 int NEMPTY, I, LENI, L, LEND, K, KLAST, ILAST, LPRINT;
522
523 NEMPTY = 0;
524 for(I = 1; I <= N; I++) {
525 LENI = LEN[I];
526 if(LENI>0) {
527 L = (LOC[I]+LENI)-1;
528 LEN[I] = IND[L];
529 IND[L] = -(N+I);
530 }
531 else if(LENI==0)
532 NEMPTY++;
533 }
534 K = 0;
535 /* Previous k */
536 KLAST = 0;
537 /* Last entry moved. */
538 ILAST = 0;
539 LEND = *LTOP;
540 for(L = 1; L <= LEND; L++) {
541 I = IND[L];
542 if(I>0) {
543 K++;
544 IND[K] = I;
545 if(REALS)
546 LUSOL->a[K] = LUSOL->a[L];
547 }
548 else if(I<-N) {
549 /* This is the end of entry i. */
550 I = -(N+I);
551 ILAST = I;
552 K++;
553 IND[K] = LEN[I];
554 if(REALS)
555 LUSOL->a[K] = LUSOL->a[L];
556 LOC[I] = KLAST+1;
557 LEN[I] = K-KLAST;
558 KLAST = K;
559 }
560 }
561 /* Move any empty items to the end, adding 1 free entry for each. */
562 if(NEMPTY>0) {
563 for(I = 1; I <= N; I++) {
564 if(LEN[I]==0) {
565 K++;
566 LOC[I] = K;
567 IND[K] = 0;
568 ILAST = I;
569 }
570 }
571 }
572 LPRINT = LUSOL->luparm[LUSOL_IP_PRINTLEVEL];
573 if(LPRINT>=LUSOL_MSG_PIVOT)
574 LUSOL_report(LUSOL, 0, "lu1rec. File compressed from %d to %d\n",
575 *LTOP,K,REALS,NEMPTY);
576 /* ncp */
577 LUSOL->luparm[LUSOL_IP_COMPRESSIONS_LU]++;
578 /* Return ilast in ind(ltop + 1). */
579 *LTOP = K;
580 IND[(*LTOP)+1] = ILAST;
581 }
582
583 /* ==================================================================
584 lu1slk sets w(j) > 0 if column j is a unit vector.
585 ------------------------------------------------------------------
586 21 Nov 2000: First version. lu1fad needs it for TCP.
587 Note that w(*) is nominally an integer array,
588 but the only spare space is the double array w(*).
589 ================================================================== */
590 void LU1SLK(LUSOLrec *LUSOL)
591 {
592 int J, LC1, LQ, LQ1, LQ2;
593
594 for(J = 1; J <= LUSOL->n; J++) {
595 LUSOL->w[J] = 0;
596 }
597 LQ1 = (LUSOL->iqloc ? LUSOL->iqloc[1] : LUSOL->n+1);
598 /* LQ1 = LUSOL->iqloc[1]; This is the original version; correction above by Yin Zhang */
599 LQ2 = LUSOL->n;
600 if(LUSOL->m>1)
601 LQ2 = LUSOL->iqloc[2]-1;
602 for(LQ = LQ1; LQ <= LQ2; LQ++) {
603 J = LUSOL->iq[LQ];
604 LC1 = LUSOL->locc[J];
605 if(fabs(LUSOL->a[LC1])==1) {
606 LUSOL->w[J] = 1;
607 }
608 }
609 }
610
611 /* ==================================================================
612 lu1gau does most of the work for each step of Gaussian elimination.
613 A multiple of the pivot column is added to each other column j
614 in the pivot row. The column list is fully updated.
615 The row list is updated if there is room, but some fill-ins may
616 remain, as indicated by ifill and jfill.
617 ------------------------------------------------------------------
618 Input:
619 ilast is the row at the end of the row list.
620 jlast is the column at the end of the column list.
621 lfirst is the first column to be processed.
622 lu + 1 is the corresponding element of U in au(*).
623 nfill keeps track of pending fill-in.
624 a(*) contains the nonzeros for each column j.
625 indc(*) contains the row indices for each column j.
626 al(*) contains the new column of L. A multiple of it is
627 used to modify each column.
628 mark(*) has been set to -1, -2, -3, ... in the rows
629 corresponding to nonzero 1, 2, 3, ... of the col of L.
630 au(*) contains the new row of U. Each nonzero gives the
631 required multiple of the column of L.
632
633 Workspace:
634 markl(*) marks the nonzeros of L actually used.
635 (A different mark, namely j, is used for each column.)
636
637 Output:
638 ilast New last row in the row list.
639 jlast New last column in the column list.
640 lfirst = 0 if all columns were completed,
641 > 0 otherwise.
642 lu returns the position of the last nonzero of U
643 actually used, in case we come back in again.
644 nfill keeps track of the total extra space needed in the
645 row file.
646 ifill(ll) counts pending fill-in for rows involved in the new
647 column of L.
648 jfill(lu) marks the first pending fill-in stored in columns
649 involved in the new row of U.
650 ------------------------------------------------------------------
651 16 Apr 1989: First version of lu1gau.
652 23 Apr 1989: lfirst, lu, nfill are now input and output
653 to allow re-entry if elimination is interrupted.
654 23 Mar 2001: Introduced ilast, jlast.
655 27 Mar 2001: Allow fill-in "in situ" if there is already room
656 up to but NOT INCLUDING the end of the
657 row or column file.
658 Seems safe way to avoid overwriting empty rows/cols
659 at the end. (May not be needed though, now that we
660 have ilast and jlast.)
661 ================================================================== */
662 void LU1GAU(LUSOLrec *LUSOL, int MELIM, int NSPARE,
663 REAL SMALL, int LPIVC1, int LPIVC2, int *LFIRST, int LPIVR2,
664 int LFREE, int MINFRE, int ILAST, int *JLAST, int *LROW, int *LCOL,
665 int *LU, int *NFILL,
666 int MARK[], REAL AL[], int MARKL[], REAL AU[], int IFILL[], int JFILL[])
667 {
668 MYBOOL ATEND;
669 int LR, J, LENJ, NFREE, LC1, LC2, NDONE, NDROP, L, I, LL, K,
670 LR1, LAST, LREP, L1, L2, LC, LENI;
671 register REAL UJ;
672 REAL AIJ;
673
674 for(LR = *LFIRST; LR <= LPIVR2; LR++) {
675 J = LUSOL->indr[LR];
676 LENJ = LUSOL->lenc[J];
677 NFREE = LFREE - *LCOL;
678 if(NFREE<MINFRE)
679 goto x900;
680 /* ---------------------------------------------------------------
681 Inner loop to modify existing nonzeros in column j.
682 Loop 440 performs most of the arithmetic involved in the
683 whole LU factorization.
684 ndone counts how many multipliers were used.
685 ndrop counts how many modified nonzeros are negligibly small.
686 --------------------------------------------------------------- */
687 (*LU)++;
688 UJ = AU[*LU];
689 LC1 = LUSOL->locc[J];
690 LC2 = (LC1+LENJ)-1;
691 ATEND = (MYBOOL) (J==*JLAST);
692 NDONE = 0;
693 if(LENJ==0)
694 goto x500;
695 NDROP = 0;
696 for(L = LC1; L <= LC2; L++) {
697 I = LUSOL->indc[L];
698 LL = -MARK[I];
699 if(LL>0) {
700 NDONE++;
701 MARKL[LL] = J;
702 LUSOL->a[L] += AL[LL]*UJ;
703 if(fabs(LUSOL->a[L])<=SMALL) {
704 NDROP++;
705 }
706 }
707 }
708 /* ---------------------------------------------------------------
709 Remove any negligible modified nonzeros from both
710 the column file and the row file.
711 --------------------------------------------------------------- */
712 if(NDROP==0)
713 goto x500;
714 K = LC1;
715 for(L = LC1; L <= LC2; L++) {
716 I = LUSOL->indc[L];
717 if(fabs(LUSOL->a[L])<=SMALL)
718 goto x460;
719 LUSOL->a[K] = LUSOL->a[L];
720 LUSOL->indc[K] = I;
721 K++;
722 continue;
723 /* Delete the nonzero from the row file. */
724 x460:
725 LENJ--;
726 LUSOL->lenr[I]--;
727 LR1 = LUSOL->locr[I];
728 LAST = LR1+LUSOL->lenr[I];
729 for(LREP = LR1; LREP <= LAST; LREP++) {
730 if(LUSOL->indr[LREP]==J)
731 break;
732 }
733 LUSOL->indr[LREP] = LUSOL->indr[LAST];
734 LUSOL->indr[LAST] = 0;
735 if(I==ILAST)
736 (*LROW)--;
737 }
738 /* Free the deleted elements from the column file. */
739 #ifdef LUSOLFastClear
740 MEMCLEAR(LUSOL->indc+K, LC2-K+1);
741 #else
742 for(L = K; L <= LC2; L++)
743 LUSOL->indc[L] = ZERO;
744 #endif
745 if(ATEND)
746 *LCOL = K-1;
747 /* ---------------------------------------------------------------
748 Deal with the fill-in in column j.
749 --------------------------------------------------------------- */
750 x500:
751 if(NDONE==MELIM)
752 goto x590;
753 /* See if column j already has room for the fill-in. */
754 if(ATEND)
755 goto x540;
756 LAST = (LC1+LENJ)-1;
757 L1 = LAST+1;
758 L2 = (LAST+MELIM)-NDONE;
759 /* 27 Mar 2001: Be sure it's not at or past end of the col file. */
760 if(L2>=*LCOL)
761 goto x520;
762 for(L = L1; L <= L2; L++) {
763 if(LUSOL->indc[L]!=0)
764 goto x520;
765 }
766 goto x540;
767 /* We must move column j to the end of the column file.
768 First, leave some spare room at the end of the
769 current last column. */
770 x520:
771 #if 1
772 L1 = (*LCOL)+1;
773 L2 = (*LCOL)+NSPARE;
774 *LCOL = L2;
775 for(L = L1; L <= L2; L++) {
776 #else
777 for(L = (*LCOL)+1; L <= (*LCOL)+NSPARE; L++) {
778 *LCOL = L; /* ****** ERROR ???? */
779 #endif
780 /* Spare space is free. */
781 LUSOL->indc[L] = 0;
782 }
783 ATEND = TRUE;
784 *JLAST = J;
785 L1 = LC1;
786 L2 = *LCOL;
787 LC1 = L2+1;
788 LUSOL->locc[J] = LC1;
789 for(L = L1; L <= LAST; L++) {
790 L2++;
791 LUSOL->a[L2] = LUSOL->a[L];
792 LUSOL->indc[L2] = LUSOL->indc[L];
793 /* Free space. */
794 LUSOL->indc[L] = 0;
795 }
796 *LCOL = L2;
797 /* ---------------------------------------------------------------
798 Inner loop for the fill-in in column j.
799 This is usually not very expensive.
800 --------------------------------------------------------------- */
801 x540:
802 LAST = (LC1+LENJ)-1;
803 LL = 0;
804 for(LC = LPIVC1; LC <= LPIVC2; LC++) {
805 LL++;
806 if(MARKL[LL]==J)
807 continue;
808 AIJ = AL[LL]*UJ;
809 if(fabs(AIJ)<=SMALL)
810 continue;
811 LENJ++;
812 LAST++;
813 LUSOL->a[LAST] = AIJ;
814 I = LUSOL->indc[LC];
815 LUSOL->indc[LAST] = I;
816 LENI = LUSOL->lenr[I];
817 /* Add 1 fill-in to row i if there is already room.
818 27 Mar 2001: Be sure it's not at or past the }
819 of the row file. */
820 L = LUSOL->locr[I]+LENI;
821 if(L>=*LROW)
822 goto x550;
823 if(LUSOL->indr[L]>0)
824 goto x550;
825 LUSOL->indr[L] = J;
826 LUSOL->lenr[I] = LENI+1;
827 continue;
828 /* Row i does not have room for the fill-in.
829 Increment ifill(ll) to count how often this has
830 happened to row i. Also, add m to the row index
831 indc(last) in column j to mark it as a fill-in that is
832 still pending.
833 If this is the first pending fill-in for row i,
834 nfill includes the current length of row i
835 (since the whole row has to be moved later).
836 If this is the first pending fill-in for column j,
837 jfill(lu) records the current length of column j
838 (to shorten the search for pending fill-ins later). */
839 x550:
840 if(IFILL[LL]==0)
841 (*NFILL) += LENI+NSPARE;
842 if(JFILL[*LU]==0)
843 JFILL[*LU] = LENJ;
844 (*NFILL)++;
845 IFILL[LL]++;
846 LUSOL->indc[LAST] = LUSOL->m+I;
847 }
848 if(ATEND)
849 *LCOL = LAST;
850 /* End loop for column j. Store its final length. */
851 x590:
852 LUSOL->lenc[J] = LENJ;
853 }
854 /* Successful completion. */
855 *LFIRST = 0;
856 return;
857 /* Interruption. We have to come back in after the
858 column file is compressed. Give lfirst a new value.
859 lu and nfill will retain their current values. */
860 x900:
861 *LFIRST = LR;
862 }
863
864 /* ==================================================================
865 lu1mar uses a Markowitz criterion to select a pivot element
866 for the next stage of a sparse LU factorization,
867 subject to a Threshold Partial Pivoting stability criterion (TPP)
868 that bounds the elements of L.
869 ------------------------------------------------------------------
870 gamma is "gamma" in the tie-breaking rule TB4 in the LUSOL paper.
871 ------------------------------------------------------------------
872 Search cols of length nz = 1, then rows of length nz = 1,
873 then cols of length nz = 2, then rows of length nz = 2, etc.
874 ------------------------------------------------------------------
875 00 Jan 1986 Version documented in LUSOL paper:
876 Gill, Murray, Saunders and Wright (1987),
877 Maintaining LU factors of a general sparse matrix,
878 Linear algebra and its applications 88/89, 239-270.
879 02 Feb 1989 Following Suhl and Aittoniemi (1987), the largest
880 element in each column is now kept at the start of
881 the column, i.e. in position locc(j) of a and indc.
882 This should speed up the Markowitz searches.
883 26 Apr 1989 Both columns and rows searched during spars1 phase.
884 Only columns searched during spars2 phase.
885 maxtie replaced by maxcol and maxrow.
886 05 Nov 1993 Initializing "mbest = m * n" wasn't big enough when
887 m = 10, n = 3, and last column had 7 nonzeros.
888 09 Feb 1994 Realised that "mbest = maxmn * maxmn" might overflow.
889 Changed to "mbest = maxmn * 1000".
890 27 Apr 2000 On large example from Todd Munson,
891 that allowed "if (mbest .le. nz1**2) go to 900"
892 to exit before any pivot had been found.
893 Introduced kbest = mbest / nz1.
894 Most pivots can be rejected with no integer multiply.
895 TRUE merit is evaluated only if it's as good as the
896 best so far (or better). There should be no danger
897 of integer overflow unless A is incredibly
898 large and dense.
899 10 Sep 2000 TCP, aijtol added for Threshold Complete Pivoting.
900 ================================================================== */
901 void LU1MAR(LUSOLrec *LUSOL, int MAXMN, MYBOOL TCP, REAL AIJTOL, REAL LTOL,
902 int MAXCOL, int MAXROW, int *IBEST, int *JBEST, int *MBEST)
903 {
904 int KBEST, NCOL, NROW, NZ1, NZ, LQ1, LQ2, LQ, J, LC1, LC2, LC, I, LEN1, MERIT, LP1,
905 LP2, LP, LR1, LR2, LR;
906 REAL ABEST, LBEST, AMAX, AIJ, CMAX;
907
908 ABEST = ZERO;
909 LBEST = ZERO;
910 *IBEST = 0;
911 *MBEST = -1;
912 KBEST = MAXMN+1;
913 NCOL = 0;
914 NROW = 0;
915 NZ1 = 0;
916 for(NZ = 1; NZ <= MAXMN; NZ++) {
917 /* nz1 = nz - 1
918 if (mbest .le. nz1**2) go to 900 */
919 if(KBEST<=NZ1)
920 goto x900;
921 if(*IBEST>0) {
922 if(NCOL>=MAXCOL)
923 goto x200;
924 }
925 if(NZ>LUSOL->m)
926 goto x200;
927 /* ---------------------------------------------------------------
928 Search the set of columns of length nz.
929 --------------------------------------------------------------- */
930 LQ1 = LUSOL->iqloc[NZ];
931 LQ2 = LUSOL->n;
932 if(NZ<LUSOL->m)
933 LQ2 = LUSOL->iqloc[NZ+1]-1;
934 for(LQ = LQ1; LQ <= LQ2; LQ++) {
935 NCOL = NCOL+1;
936 J = LUSOL->iq[LQ];
937 LC1 = LUSOL->locc[J];
938 LC2 = LC1+NZ1;
939 AMAX = fabs(LUSOL->a[LC1]);
940 /* Test all aijs in this column.
941 amax is the largest element (the first in the column).
942 cmax is the largest multiplier if aij becomes pivot. */
943 if(TCP) {
944 /* Nothing in whole column */
945 if(AMAX<AIJTOL)
946 continue;
947 }
948 for(LC = LC1; LC <= LC2; LC++) {
949 I = LUSOL->indc[LC];
950 LEN1 = LUSOL->lenr[I]-1;
951 /* merit = nz1 * len1
952 if (merit > mbest) continue; */
953 if(LEN1>KBEST)
954 continue;
955 /* aij has a promising merit.
956 Apply the stability test.
957 We require aij to be sufficiently large compared to
958 all other nonzeros in column j. This is equivalent
959 to requiring cmax to be bounded by Ltol. */
960 if(LC==LC1) {
961 /* This is the maximum element, amax.
962 Find the biggest element in the rest of the column
963 and hence get cmax. We know cmax .le. 1, but
964 we still want it exactly in order to break ties.
965 27 Apr 2002: Settle for cmax = 1. */
966 AIJ = AMAX;
967 CMAX = ONE;
968 /* cmax = zero
969 for (l = lc1 + 1; l <= lc2; l++)
970 cmax = max( cmax, abs( a(l) ) );
971 cmax = cmax / amax; */
972 }
973 else {
974 /* aij is not the biggest element, so cmax .ge. 1.
975 Bail out if cmax will be too big. */
976 AIJ = fabs(LUSOL->a[LC]);
977 /* Absolute test for Complete Pivoting */
978 if(TCP) {
979 if(AIJ<AIJTOL)
980 continue;
981 /* TPP */
982 }
983 else {
984 if(AIJ*LTOL<AMAX)
985 continue;
986 }
987 CMAX = AMAX/AIJ;
988 }
989 /* aij is big enough. Its maximum multiplier is cmax. */
990 MERIT = NZ1*LEN1;
991 if(MERIT==*MBEST) {
992 /* Break ties.
993 (Initializing mbest < 0 prevents getting here if
994 nothing has been found yet.)
995 In this version we minimize cmax
996 but if it is already small we maximize the pivot. */
997 if(LBEST<=LUSOL->parmlu[LUSOL_RP_GAMMA] &&
998 CMAX<=LUSOL->parmlu[LUSOL_RP_GAMMA]) {
999 if(ABEST>=AIJ)
1000 continue;
1001 }
1002 else {
1003 if(LBEST<=CMAX)
1004 continue;
1005 }
1006 }
1007 /* aij is the best pivot so far. */
1008 *IBEST = I;
1009 *JBEST = J;
1010 KBEST = LEN1;
1011 *MBEST = MERIT;
1012 ABEST = AIJ;
1013 LBEST = CMAX;
1014 if(NZ==1)
1015 goto x900;
1016 }
1017 /* Finished with that column. */
1018 if(*IBEST>0) {
1019 if(NCOL>=MAXCOL)
1020 goto x200;
1021 }
1022 }
1023 /* ---------------------------------------------------------------
1024 Search the set of rows of length nz.
1025 --------------------------------------------------------------- */
1026 x200:
1027 /* if (mbest .le. nz*nz1) go to 900 */
1028 if(KBEST<=NZ)
1029 goto x900;
1030 if(*IBEST>0) {
1031 if(NROW>=MAXROW)
1032 goto x290;
1033 }
1034 if(NZ>LUSOL->n)
1035 goto x290;
1036 LP1 = LUSOL->iploc[NZ];
1037 LP2 = LUSOL->m;
1038 if(NZ<LUSOL->n)
1039 LP2 = LUSOL->iploc[NZ+1]-1;
1040 for(LP = LP1; LP <= LP2; LP++) {
1041 NROW++;
1042 I = LUSOL->ip[LP];
1043 LR1 = LUSOL->locr[I];
1044 LR2 = LR1+NZ1;
1045 for(LR = LR1; LR <= LR2; LR++) {
1046 J = LUSOL->indr[LR];
1047 LEN1 = LUSOL->lenc[J]-1;
1048 /* merit = nz1 * len1
1049 if (merit .gt. mbest) continue */
1050 if(LEN1>KBEST)
1051 continue;
1052 /* aij has a promising merit.
1053 Find where aij is in column j. */
1054 LC1 = LUSOL->locc[J];
1055 LC2 = LC1+LEN1;
1056 AMAX = fabs(LUSOL->a[LC1]);
1057 for(LC = LC1; LC <= LC2; LC++) {
1058 if(LUSOL->indc[LC]==I)
1059 break;
1060 }
1061 /* Apply the same stability test as above. */
1062 AIJ = fabs(LUSOL->a[LC]);
1063 /* Absolute test for Complete Pivoting */
1064 if(TCP) {
1065 if(AIJ<AIJTOL)
1066 continue;
1067 }
1068 if(LC==LC1) {
1069 /* This is the maximum element, amax.
1070 Find the biggest element in the rest of the column
1071 and hence get cmax. We know cmax .le. 1, but
1072 we still want it exactly in order to break ties.
1073 27 Apr 2002: Settle for cmax = 1. */
1074 CMAX = ONE;
1075 /* cmax = zero
1076 for(l = lc1 + 1; l <= lc2; l++)
1077 cmax = max( cmax, fabs( a(l) ) )
1078 cmax = cmax / amax */
1079 }
1080 else {
1081 /* aij is not the biggest element, so cmax .ge. 1.
1082 Bail out if cmax will be too big. */
1083 if(TCP) {
1084 /* relax */
1085 }
1086 else {
1087 if(AIJ*LTOL<AMAX)
1088 continue;
1089 }
1090 CMAX = AMAX/AIJ;
1091 }
1092 /* aij is big enough. Its maximum multiplier is cmax. */
1093 MERIT = NZ1*LEN1;
1094 if(MERIT==*MBEST) {
1095 /* Break ties as before.
1096 (Initializing mbest < 0 prevents getting here if
1097 nothing has been found yet.) */
1098 if(LBEST<=LUSOL->parmlu[LUSOL_RP_GAMMA] &&
1099 CMAX<=LUSOL->parmlu[LUSOL_RP_GAMMA]) {
1100 if(ABEST>=AIJ)
1101 continue;
1102 }
1103 else {
1104 if(LBEST<=CMAX)
1105 continue;
1106 }
1107 }
1108 /* aij is the best pivot so far. */
1109 *IBEST = I;
1110 *JBEST = J;
1111 *MBEST = MERIT;
1112 KBEST = LEN1;
1113 ABEST = AIJ;
1114 LBEST = CMAX;
1115 if(NZ==1)
1116 goto x900;
1117 }
1118 /* Finished with that row. */
1119 if(*IBEST>0) {
1120 if(NROW>=MAXROW)
1121 goto x290;
1122 }
1123 }
1124 /* See if it's time to quit. */
1125 x290:
1126 if(*IBEST>0) {
1127 if(NROW>=MAXROW && NCOL>=MAXCOL)
1128 goto x900;
1129 }
1130 /* Press on with next nz. */
1131 NZ1 = NZ;
1132 if(*IBEST>0)
1133 KBEST = *MBEST/NZ1;
1134 }
1135 x900:
1136 ;
1137 }
1138
1139 /* ==================================================================
1140 lu1mCP uses a Markowitz criterion to select a pivot element
1141 for the next stage of a sparse LU factorization,
1142 subject to a Threshold Complete Pivoting stability criterion (TCP)
1143 that bounds the elements of L and U.
1144 ------------------------------------------------------------------
1145 gamma is "gamma" in the tie-breaking rule TB4 in the LUSOL paper.
1146 ------------------------------------------------------------------
1147 09 May 2002: First version of lu1mCP.
1148 It searches columns only, using the heap that
1149 holds the largest element in each column.
1150 09 May 2002: Current version of lu1mCP.
1151 ================================================================== */
1152 void LU1MCP(LUSOLrec *LUSOL, REAL AIJTOL, int *IBEST, int *JBEST, int *MBEST,
1153 int HLEN, REAL HA[], int HJ[])
1154 {
1155 int J, KHEAP, LC, LC1, LC2, LENJ, MAXCOL, NCOL, NZ1, I, LEN1, MERIT;
1156 REAL ABEST, AIJ, AMAX, CMAX, LBEST;
1157
1158 /* ------------------------------------------------------------------
1159 Search up to maxcol columns stored at the top of the heap.
1160 The very top column helps initialize mbest.
1161 ------------------------------------------------------------------ */
1162 ABEST = ZERO;
1163 LBEST = ZERO;
1164 *IBEST = 0;
1165 /* Column at the top of the heap */
1166 *JBEST = HJ[1];
1167 LENJ = LUSOL->lenc[*JBEST];
1168 /* Bigger than any possible merit */
1169 *MBEST = LENJ*HLEN;
1170 /* ??? Big question */
1171 MAXCOL = 40;
1172 /* No. of columns searched */
1173 NCOL = 0;
1174 for(KHEAP = 1; KHEAP <= HLEN; KHEAP++) {
1175 AMAX = HA[KHEAP];
1176 if(AMAX<AIJTOL)
1177 continue;
1178 NCOL++;
1179 J = HJ[KHEAP];
1180 /* ---------------------------------------------------------------
1181 This column has at least one entry big enough (the top one).
1182 Search the column for other possibilities.
1183 --------------------------------------------------------------- */
1184 LENJ = LUSOL->lenc[J];
1185 NZ1 = LENJ-1;
1186 LC1 = LUSOL->locc[J];
1187 LC2 = LC1+NZ1;
1188 /* -- amax = abs( a(lc1) )
1189 Test all aijs in this column.
1190 amax is the largest element (the first in the column).
1191 cmax is the largest multiplier if aij becomes pivot. */
1192 for(LC = LC1; LC <= LC2; LC++) {
1193 I = LUSOL->indc[LC];
1194 LEN1 = LUSOL->lenr[I]-1;
1195 MERIT = NZ1*LEN1;
1196 if(MERIT>*MBEST)
1197 continue;
1198 /* aij has a promising merit. */
1199 if(LC==LC1) {
1200 /* This is the maximum element, amax.
1201 Find the biggest element in the rest of the column
1202 and hence get cmax. We know cmax .le. 1, but
1203 we still want it exactly in order to break ties.
1204 27 Apr 2002: Settle for cmax = 1. */
1205 AIJ = AMAX;
1206 CMAX = ONE;
1207 /* cmax = ZERO;
1208 for(l = lc1 + 1; l <= lc2; l++)
1209 cmax = max( cmax, abs( a(l) ) )
1210 cmax = cmax / amax; */
1211 }
1212 else {
1213 /* aij is not the biggest element, so cmax .ge. 1.
1214 Bail out if cmax will be too big. */
1215 AIJ = fabs(LUSOL->a[LC]);
1216 if(AIJ<AIJTOL)
1217 continue;
1218 CMAX = AMAX/AIJ;
1219 }
1220 /* aij is big enough. Its maximum multiplier is cmax. */
1221 if(MERIT==*MBEST) {
1222 /* Break ties.
1223 (Initializing mbest "too big" prevents getting here if
1224 nothing has been found yet.)
1225 In this version we minimize cmax
1226 but if it is already small we maximize the pivot. */
1227 if(LBEST<=LUSOL->parmlu[LUSOL_RP_GAMMA] &&
1228 CMAX<=LUSOL->parmlu[LUSOL_RP_GAMMA]) {
1229 if(ABEST>=AIJ)
1230 continue;
1231 }
1232 else {
1233 if(LBEST<=CMAX)
1234 continue;
1235 }
1236 }
1237 /* aij is the best pivot so far. */
1238 *IBEST = I;
1239 *JBEST = J;
1240 *MBEST = MERIT;
1241 ABEST = AIJ;
1242 LBEST = CMAX;
1243 /* Col or row of length 1 */
1244 if(MERIT==0)
1245 goto x900;
1246 }
1247 if(NCOL>=MAXCOL)
1248 goto x900;
1249 }
1250 x900:
1251 ;
1252 }
1253
1254 /* ==================================================================
1255 lu1mRP uses a Markowitz criterion to select a pivot element
1256 for the next stage of a sparse LU factorization,
1257 subject to a Threshold Rook Pivoting stability criterion (TRP)
1258 that bounds the elements of L and U.
1259 ------------------------------------------------------------------
1260 11 Jun 2002: First version of lu1mRP derived from lu1mar.
1261 11 Jun 2002: Current version of lu1mRP.
1262 ================================================================== */
1263 void LU1MRP(LUSOLrec *LUSOL, int MAXMN, REAL LTOL, int MAXCOL, int MAXROW,
1264 int *IBEST, int *JBEST, int *MBEST, REAL AMAXR[])
1265 {
1266 int I, J, KBEST, LC, LC1, LC2, LEN1, LP, LP1, LP2, LQ, LQ1,
1267 LQ2, LR, LR1, LR2, MERIT, NCOL, NROW, NZ, NZ1;
1268 REAL ABEST, AIJ, AMAX, ATOLI, ATOLJ;
1269
1270 /* ------------------------------------------------------------------
1271 Search cols of length nz = 1, then rows of length nz = 1,
1272 then cols of length nz = 2, then rows of length nz = 2, etc.
1273 ------------------------------------------------------------------ */
1274 ABEST = ZERO;
1275 *IBEST = 0;
1276 KBEST = MAXMN+1;
1277 *MBEST = -1;
1278 NCOL = 0;
1279 NROW = 0;
1280 NZ1 = 0;
1281 for(NZ = 1; NZ <= MAXMN; NZ++) {
1282 /* nz1 = nz - 1
1283 if (mbest .le. nz1**2) go to 900 */
1284 if(KBEST<=NZ1)
1285 goto x900;
1286 if(*IBEST>0) {
1287 if(NCOL>=MAXCOL)
1288 goto x200;
1289 }
1290 if(NZ>LUSOL->m)
1291 goto x200;
1292 /* ---------------------------------------------------------------
1293 Search the set of columns of length nz.
1294 --------------------------------------------------------------- */
1295 LQ1 = LUSOL->iqloc[NZ];
1296 LQ2 = LUSOL->n;
1297 if(NZ<LUSOL->m)
1298 LQ2 = LUSOL->iqloc[NZ+1]-1;
1299 for(LQ = LQ1; LQ <= LQ2; LQ++) {
1300 NCOL = NCOL+1;
1301 J = LUSOL->iq[LQ];
1302 LC1 = LUSOL->locc[J];
1303 LC2 = LC1+NZ1;
1304 AMAX = fabs(LUSOL->a[LC1]);
1305 /* Min size of pivots in col j */
1306 ATOLJ = AMAX/LTOL;
1307 /* Test all aijs in this column. */
1308 for(LC = LC1; LC <= LC2; LC++) {
1309 I = LUSOL->indc[LC];
1310 LEN1 = LUSOL->lenr[I]-1;
1311 /* merit = nz1 * len1
1312 if (merit .gt. mbest) continue; */
1313 if(LEN1>KBEST)
1314 continue;
1315 /* aij has a promising merit.
1316 Apply the Threshold Rook Pivoting stability test.
1317 First we require aij to be sufficiently large
1318 compared to other nonzeros in column j.
1319 Then we require aij to be sufficiently large
1320 compared to other nonzeros in row i. */
1321 AIJ = fabs(LUSOL->a[LC]);
1322 if(AIJ<ATOLJ)
1323 continue;
1324 if(AIJ*LTOL<AMAXR[I])
1325 continue;
1326 /* aij is big enough. */
1327 MERIT = NZ1*LEN1;
1328 if(MERIT==*MBEST) {
1329 /* Break ties.
1330 (Initializing mbest < 0 prevents getting here if
1331 nothing has been found yet.) */
1332 if(ABEST>=AIJ)
1333 continue;
1334 }
1335 /* aij is the best pivot so far. */
1336 *IBEST = I;
1337 *JBEST = J;
1338 KBEST = LEN1;
1339 *MBEST = MERIT;
1340 ABEST = AIJ;
1341 if(NZ==1)
1342 goto x900;
1343 }
1344 /* Finished with that column. */
1345 if(*IBEST>0) {
1346 if(NCOL>=MAXCOL)
1347 goto x200;
1348 }
1349 }
1350 /* ---------------------------------------------------------------
1351 Search the set of rows of length nz.
1352 --------------------------------------------------------------- */
1353 x200:
1354 /* if (mbest .le. nz*nz1) go to 900 */
1355 if(KBEST<=NZ)
1356 goto x900;
1357 if(*IBEST>0) {
1358 if(NROW>=MAXROW)
1359 goto x290;
1360 }
1361 if(NZ>LUSOL->n)
1362 goto x290;
1363 LP1 = LUSOL->iploc[NZ];
1364 LP2 = LUSOL->m;
1365 if(NZ<LUSOL->n)
1366 LP2 = LUSOL->iploc[NZ+1]-1;
1367 for(LP = LP1; LP <= LP2; LP++) {
1368 NROW = NROW+1;
1369 I = LUSOL->ip[LP];
1370 LR1 = LUSOL->locr[I];
1371 LR2 = LR1+NZ1;
1372 /* Min size of pivots in row i */
1373 ATOLI = AMAXR[I]/LTOL;
1374 for(LR = LR1; LR <= LR2; LR++) {
1375 J = LUSOL->indr[LR];
1376 LEN1 = LUSOL->lenc[J]-1;
1377 /* merit = nz1 * len1
1378 if (merit .gt. mbest) continue; */
1379 if(LEN1>KBEST)
1380 continue;
1381 /* aij has a promising merit.
1382 Find where aij is in column j. */
1383 LC1 = LUSOL->locc[J];
1384 LC2 = LC1+LEN1;
1385 AMAX = fabs(LUSOL->a[LC1]);
1386 for(LC = LC1; LC <= LC2; LC++) {
1387 if(LUSOL->indc[LC]==I)
1388 break;
1389 }
1390 /* Apply the Threshold Rook Pivoting stability test.
1391 First we require aij to be sufficiently large
1392 compared to other nonzeros in row i.
1393 Then we require aij to be sufficiently large
1394 compared to other nonzeros in column j. */
1395 AIJ = fabs(LUSOL->a[LC]);
1396 if(AIJ<ATOLI)
1397 continue;
1398 if(AIJ*LTOL<AMAX)
1399 continue;
1400 /* aij is big enough. */
1401 MERIT = NZ1*LEN1;
1402 if(MERIT==*MBEST) {
1403 /* Break ties as before.
1404 (Initializing mbest < 0 prevents getting here if
1405 nothing has been found yet.) */
1406 if(ABEST>=AIJ)
1407 continue;
1408 }
1409 /* aij is the best pivot so far. */
1410 *IBEST = I;
1411 *JBEST = J;
1412 KBEST = LEN1;
1413 *MBEST = MERIT;
1414 ABEST = AIJ;
1415 if(NZ==1)
1416 goto x900;
1417 }
1418 /* Finished with that row. */
1419 if(*IBEST>0) {
1420 if(NROW>=MAXROW)
1421 goto x290;
1422 }
1423 }
1424 /* See if it's time to quit. */
1425 x290:
1426 if(*IBEST>0) {
1427 if(NROW>=MAXROW && NCOL>=MAXCOL)
1428 goto x900;
1429 }
1430 /* Press on with next nz. */
1431 NZ1 = NZ;
1432 if(*IBEST>0)
1433 KBEST = *MBEST/NZ1;
1434 }
1435 x900:
1436 ;
1437 }
1438
1439 /* ==================================================================
1440 lu1mSP is intended for symmetric matrices that are either
1441 definite or quasi-definite.
1442 lu1mSP uses a Markowitz criterion to select a pivot element for
1443 the next stage of a sparse LU factorization of a symmetric matrix,
1444 subject to a Threshold Symmetric Pivoting stability criterion
1445 (TSP) restricted to diagonal elements to preserve symmetry.
1446 This bounds the elements of L and U and should have rank-revealing
1447 properties analogous to Threshold Rook Pivoting for unsymmetric
1448 matrices.
1449 ------------------------------------------------------------------
1450 14 Dec 2002: First version of lu1mSP derived from lu1mRP.
1451 There is no safeguard to ensure that A is symmetric.
1452 14 Dec 2002: Current version of lu1mSP.
1453 ================================================================== */
1454 void LU1MSP(LUSOLrec *LUSOL, int MAXMN, REAL LTOL, int MAXCOL,
1455 int *IBEST, int *JBEST, int *MBEST)
1456 {
1457 int I, J, KBEST, LC, LC1, LC2, LQ, LQ1, LQ2, MERIT, NCOL, NZ, NZ1;
1458 REAL ABEST, AIJ, AMAX, ATOLJ;
1459
1460 /* ------------------------------------------------------------------
1461 Search cols of length nz = 1, then cols of length nz = 2, etc.
1462 ------------------------------------------------------------------ */
1463 ABEST = ZERO;
1464 *IBEST = 0;
1465 *MBEST = -1;
1466 KBEST = MAXMN+1;
1467 NCOL = 0;
1468 NZ1 = 0;
1469 for(NZ = 1; NZ <= MAXMN; NZ++) {
1470 /* nz1 = nz - 1
1471 if (mbest .le. nz1**2) go to 900 */
1472 if(KBEST<=NZ1)
1473 goto x900;
1474 if(*IBEST>0) {
1475 if(NCOL>=MAXCOL)
1476 goto x200;
1477 }
1478 if(NZ>LUSOL->m)
1479 goto x200;
1480 /* ---------------------------------------------------------------
1481 Search the set of columns of length nz.
1482 --------------------------------------------------------------- */
1483 LQ1 = LUSOL->iqloc[NZ];
1484 LQ2 = LUSOL->n;
1485 if(NZ<LUSOL->m)
1486 LQ2 = LUSOL->iqloc[NZ+1]-1;
1487 for(LQ = LQ1; LQ <= LQ2; LQ++) {
1488 NCOL++;
1489 J = LUSOL->iq[LQ];
1490 LC1 = LUSOL->locc[J];
1491 LC2 = LC1+NZ1;
1492 AMAX = fabs(LUSOL->a[LC1]);
1493 /* Min size of pivots in col j */
1494 ATOLJ = AMAX/LTOL;
1495 /* Test all aijs in this column.
1496 Ignore everything except the diagonal. */
1497 for(LC = LC1; LC <= LC2; LC++) {
1498 I = LUSOL->indc[LC];
1499 /* Skip off-diagonals. */
1500 if(I!=J)
1501 continue;
1502 /* merit = nz1 * nz1
1503 if (merit .gt. mbest) continue; */
1504 if(NZ1>KBEST)
1505 continue;
1506 /* aij has a promising merit.
1507 Apply the Threshold Partial Pivoting stability test
1508 (which is equivalent to Threshold Rook Pivoting for
1509 symmetric matrices).
1510 We require aij to be sufficiently large
1511 compared to other nonzeros in column j. */
1512 AIJ = fabs(LUSOL->a[LC]);
1513 if(AIJ<ATOLJ)
1514 continue;
1515 /* aij is big enough. */
1516 MERIT = NZ1*NZ1;
1517 if(MERIT==*MBEST) {
1518 /* Break ties.
1519 (Initializing mbest < 0 prevents getting here if
1520 nothing has been found yet.) */
1521 if(ABEST>=AIJ)
1522 continue;
1523 }
1524 /* aij is the best pivot so far. */
1525 *IBEST = I;
1526 *JBEST = J;
1527 KBEST = NZ1;
1528 *MBEST = MERIT;
1529 ABEST = AIJ;
1530 if(NZ==1)
1531 goto x900;
1532 }
1533 /* Finished with that column. */
1534 if(*IBEST>0) {
1535 if(NCOL>=MAXCOL)
1536 goto x200;
1537 }
1538 }
1539 /* See if it's time to quit. */
1540 x200:
1541 if(*IBEST>0) {
1542 if(NCOL>=MAXCOL)
1543 goto x900;
1544 }
1545 /* Press on with next nz. */
1546 NZ1 = NZ;
1547 if(*IBEST>0)
1548 KBEST = *MBEST/NZ1;
1549 }
1550 x900:
1551 ;
1552 }
1553
1554 /* ==================================================================
1555 lu1mxc moves the largest element in each of columns iq(k1:k2)
1556 to the top of its column.
1557 If k1 > k2, nothing happens.
1558 ------------------------------------------------------------------
1559 06 May 2002: (and earlier)
1560 All columns k1:k2 must have one or more elements.
1561 07 May 2002: Allow for empty columns. The heap routines need to
1562 find 0.0 as the "largest element".
1563 29 Nov 2005: Bug fix - avoiding overwriting the next column when
1564 the current column is empty (i.e. LENJ==0)
1565 Yin Zhang <yzhang@cs.utexas.edu>
1566 ================================================================== */
1567 void LU1MXC(LUSOLrec *LUSOL, int K1, int K2, int IX[])
1568 {
1569 int I, J, K, L, LC, LENJ;
1570 REAL AMAX;
1571
1572 for(K = K1; K <= K2; K++) {
1573 J = IX[K];
1574 LC = LUSOL->locc[J];
1575 LENJ = LUSOL->lenc[J];
1576 if(LENJ==0)
1577 /* LUSOL->a[LC] = ZERO; Removal suggested by Yin Zhang to avoid overwriting next column when current is empty */
1578 ;
1579 else {
1580 L = idamaxlpsolve(LUSOL->lenc[J], LUSOL->a + LC - LUSOL_ARRAYOFFSET,1) + LC - 1;
1581 if(L>LC) {
1582 AMAX = LUSOL->a[L];
1583 LUSOL->a[L] = LUSOL->a[LC];
1584 LUSOL->a[LC] = AMAX;
1585 I = LUSOL->indc[L];
1586 LUSOL->indc[L] = LUSOL->indc[LC];
1587 LUSOL->indc[LC] = I;
1588 }
1589 }
1590 }
1591 }
1592
1593 /* ==================================================================
1594 lu1mxr finds the largest element in each of row ip(k1:k2)
1595 and stores it in Amaxr(*). The nonzeros are stored column-wise
1596 in (a,indc,lenc,locc) and their structure is row-wise
1597 in ( indr,lenr,locr).
1598 If k1 > k2, nothing happens.
1599 ------------------------------------------------------------------
1600 11 Jun 2002: First version of lu1mxr.
1601 Allow for empty columns.
1602 ================================================================== */
1603 void LU1MXR(LUSOLrec *LUSOL, int K1, int K2, int IX[], REAL AMAXR[])
1604 {
1605 #define FastMXR
1606 #ifdef FastMXR
1607 static int I, *J, *IC, K, LC, LC1, LC2, LR, LR1, LR2;
1608 static REAL AMAX;
1609 #else
1610 int I, J, K, LC, LC1, LC2, LR, LR1, LR2;
1611 REAL AMAX;
1612 #endif
1613
1614 for(K = K1; K <= K2; K++) {
1615 AMAX = ZERO;
1616 I = IX[K];
1617 /* Find largest element in row i. */
1618 LR1 = LUSOL->locr[I];
1619 LR2 = (LR1+LUSOL->lenr[I])-1;
1620 #ifdef FastMXR
1621 for(LR = LR1, J = LUSOL->indr + LR1;
1622 LR <= LR2; LR++, J++) {
1623 /* Find where aij is in column j. */
1624 LC1 = LUSOL->locc[*J];
1625 LC2 = LC1+LUSOL->lenc[*J];
1626 for(LC = LC1, IC = LUSOL->indc + LC1;
1627 LC < LC2; LC++, IC++) {
1628 if(*IC==I)
1629 break;
1630 }
1631 SETMAX(AMAX,fabs(LUSOL->a[LC]));
1632 }
1633 #else
1634 for(LR = LR1; LR <= LR2; LR++) {
1635 J = LUSOL->indr[LR];
1636 /* Find where aij is in column j. */
1637 LC1 = LUSOL->locc[J];
1638 LC2 = (LC1+LUSOL->lenc[J])-1;
1639 for(LC = LC1; LC <= LC2; LC++) {
1640 if(LUSOL->indc[LC]==I)
1641 break;
1642 }
1643 SETMAX(AMAX,fabs(LUSOL->a[LC]));
1644 }
1645 #endif
1646 AMAXR[I] = AMAX;
1647 }
1648 }
1649
1650
1651 /* ==================================================================
1652 lu1ful computes a dense (full) LU factorization of the
1653 mleft by nleft matrix that remains to be factored at the
1654 beginning of the nrowu-th pass through the main loop of lu1fad.
1655 ------------------------------------------------------------------
1656 02 May 1989: First version.
1657 05 Feb 1994: Column interchanges added to lu1DPP.
1658 08 Feb 1994: ipinv reconstructed, since lu1pq3 may alter ip.
1659 ================================================================== */
1660 void LU1FUL(LUSOLrec *LUSOL, int LEND, int LU1, MYBOOL TPP,
1661 int MLEFT, int NLEFT, int NRANK, int NROWU,
1662 int *LENL, int *LENU, int *NSING,
1663 MYBOOL KEEPLU, REAL SMALL, REAL D[], int IPVT[])
1664 {
1665 int L, I, J, IPBASE, LDBASE, LQ, LC1, LC2, LC, LD, LKK, LKN, LU, K, L1,
1666 L2, IBEST, JBEST, LA, LL, NROWD, NCOLD;
1667 REAL AI, AJ;
1668
1669 /* ------------------------------------------------------------------
1670 If lu1pq3 moved any empty rows, reset ipinv = inverse of ip.
1671 ------------------------------------------------------------------ */
1672 if(NRANK<LUSOL->m) {
1673 for(L = 1; L <= LUSOL->m; L++) {
1674 I = LUSOL->ip[L];
1675 LUSOL->ipinv[I] = L;
1676 }
1677 }
1678 /* ------------------------------------------------------------------
1679 Copy the remaining matrix into the dense matrix D.
1680 ------------------------------------------------------------------ */
1681 #ifdef LUSOLFastClear
1682 MEMCLEAR((D+1), LEND);
1683 #else
1684 /* dload(LEND, ZERO, D, 1); */
1685 for(J = 1; J <= LEND; J++)
1686 D[J] = ZERO;
1687 #endif
1688
1689 IPBASE = NROWU-1;
1690 LDBASE = 1-NROWU;
1691 for(LQ = NROWU; LQ <= LUSOL->n; LQ++) {
1692 J = LUSOL->iq[LQ];
1693 LC1 = LUSOL->locc[J];
1694 LC2 = (LC1+LUSOL->lenc[J])-1;
1695 for(LC = LC1; LC <= LC2; LC++) {
1696 I = LUSOL->indc[LC];
1697 LD = LDBASE+LUSOL->ipinv[I];
1698 D[LD] = LUSOL->a[LC];
1699 }
1700 LDBASE += MLEFT;
1701 }
1702 /* ------------------------------------------------------------------
1703 Call our favorite dense LU factorizer.
1704 ------------------------------------------------------------------ */
1705 if(TPP)
1706 LU1DPP(LUSOL, D,MLEFT,MLEFT,NLEFT,SMALL,NSING,IPVT,LUSOL->iq+NROWU-LUSOL_ARRAYOFFSET);
1707 else
1708 LU1DCP(LUSOL, D,MLEFT,MLEFT,NLEFT,SMALL,NSING,IPVT,LUSOL->iq+NROWU-LUSOL_ARRAYOFFSET);
1709
1710 /* ------------------------------------------------------------------
1711 Move D to the beginning of A,
1712 and pack L and U at the top of a, indc, indr.
1713 In the process, apply the row permutation to ip.
1714 lkk points to the diagonal of U.
1715 ------------------------------------------------------------------ */
1716 #ifdef LUSOLFastCopy
1717 MEMCOPY(LUSOL->a+1,D+1,LEND);
1718 #else
1719 dcopylpsolve(LEND,D,1,LUSOL->a,1);
1720 #endif
1721 #ifdef ClassicdiagU
1722 LUSOL->diagU = LUSOL->a + (LUSOL->lena-LUSOL->n);
1723 #endif
1724 LKK = 1;
1725 LKN = (LEND-MLEFT)+1;
1726 LU = LU1;
1727 for(K = 1; K <= MIN(MLEFT,NLEFT); K++) {
1728 L1 = IPBASE+K;
1729 L2 = IPBASE+IPVT[K];
1730 if(L1!=L2) {
1731 I = LUSOL->ip[L1];
1732 LUSOL->ip[L1] = LUSOL->ip[L2];
1733 LUSOL->ip[L2] = I;
1734 }
1735 IBEST = LUSOL->ip[L1];
1736 JBEST = LUSOL->iq[L1];
1737 if(KEEPLU) {
1738 /* ===========================================================
1739 Pack the next column of L.
1740 =========================================================== */
1741 LA = LKK;
1742 LL = LU;
1743 NROWD = 1;
1744 for(I = K+1; I <= MLEFT; I++) {
1745 LA++;
1746 AI = LUSOL->a[LA];
1747 if(fabs(AI)>SMALL) {
1748 NROWD = NROWD+1;
1749 LL--;
1750 LUSOL->a[LL] = AI;
1751 LUSOL->indc[LL] = LUSOL->ip[IPBASE+I];
1752 LUSOL->indr[LL] = IBEST;
1753 }
1754 }
1755 /* ===========================================================
1756 Pack the next row of U.
1757 We go backwards through the row of D
1758 so the diagonal ends up at the front of the row of U.
1759 Beware -- the diagonal may be zero.
1760 =========================================================== */
1761 LA = LKN+MLEFT;
1762 LU = LL;
1763 NCOLD = 0;
1764 for(J = NLEFT; J >= K; J--) {
1765 LA = LA-MLEFT;
1766 AJ = LUSOL->a[LA];
1767 if(fabs(AJ)>SMALL || J==K) {
1768 NCOLD++;
1769 LU--;
1770 LUSOL->a[LU] = AJ;
1771 LUSOL->indr[LU] = LUSOL->iq[IPBASE+J];
1772 }
1773 }
1774 LUSOL->lenr[IBEST] = -NCOLD;
1775 LUSOL->lenc[JBEST] = -NROWD;
1776 *LENL = ((*LENL)+NROWD)-1;
1777 *LENU = (*LENU)+NCOLD;
1778 LKN++;
1779 }
1780 else {
1781 /* ===========================================================
1782 Store just the diagonal of U, in natural order.
1783 =========================================================== */
1784 LUSOL->diagU[JBEST] = LUSOL->a[LKK];
1785 }
1786 LKK += MLEFT+1;
1787 }
1788 }
1789
1790
1791 /* ==================================================================
1792 lu1or1 organizes the elements of an m by n matrix A as
1793 follows. On entry, the parallel arrays a, indc, indr,
1794 contain nelem entries of the form aij, i, j,
1795 in any order. nelem must be positive.
1796 Entries not larger than the input parameter small are treated as
1797 zero and removed from a, indc, indr. The remaining entries are
1798 defined to be nonzero. numnz returns the number of such nonzeros
1799 and Amax returns the magnitude of the largest nonzero.
1800 The arrays lenc, lenr return the number of nonzeros in each
1801 column and row of A.
1802 inform = 0 on exit, except inform = 1 if any of the indices in
1803 indc, indr imply that the element aij lies outside the m by n
1804 dimensions of A.
1805 ------------------------------------------------------------------
1806 xx Feb 1985: Original version.
1807 17 Oct 2000: a, indc, indr now have size lena to allow nelem = 0.
1808 ================================================================== */
1809 void LU1OR1(LUSOLrec *LUSOL, REAL SMALL,
1810 REAL *AMAX, int *NUMNZ, int *LERR, int *INFORM)
1811 {
1812 int I, J, L, LDUMMY;
1813
1814 #ifdef LUSOLFastClear
1815 MEMCLEAR((LUSOL->lenr+1), LUSOL->m);
1816 MEMCLEAR((LUSOL->lenc+1), LUSOL->n);
1817 #else
1818 for(I = 1; I <= LUSOL->m; I++)
1819 LUSOL->lenr[I] = ZERO;
1820 for(I = 1; I <= LUSOL->n; I++)
1821 LUSOL->lenc[I] = ZERO;
1822 #endif
1823
1824 *AMAX = 0;
1825 *NUMNZ = LUSOL->nelem;
1826 L = LUSOL->nelem+1;
1827 for(LDUMMY = 1; LDUMMY <= LUSOL->nelem; LDUMMY++) {
1828 L--;
1829 if(fabs(LUSOL->a[L])>SMALL) {
1830 I = LUSOL->indc[L];
1831 J = LUSOL->indr[L];
1832 SETMAX(*AMAX,fabs(LUSOL->a[L]));
1833 if(I<1 || I>LUSOL->m)
1834 goto x910;
1835 if(J<1 || J>LUSOL->n)
1836 goto x910;
1837 LUSOL->lenr[I]++;
1838 LUSOL->lenc[J]++;
1839 }
1840 else {
1841 /* Replace a negligible element by last element. Since
1842 we are going backwards, we know the last element is ok. */
1843 LUSOL->a[L] = LUSOL->a[*NUMNZ];
1844 LUSOL->indc[L] = LUSOL->indc[*NUMNZ];
1845 LUSOL->indr[L] = LUSOL->indr[*NUMNZ];
1846 (*NUMNZ)--;
1847 }
1848 }
1849 *LERR = 0;
1850 *INFORM = LUSOL_INFORM_LUSUCCESS;
1851 return;
1852
1853 x910:
1854 *LERR = L;
1855 *INFORM = LUSOL_INFORM_LUSINGULAR;
1856 }
1857
1858 /* ==================================================================
1859 lu1or2 sorts a list of matrix elements a(i,j) into column
1860 order, given numa entries a(i,j), i, j in the parallel
1861 arrays a, inum, jnum respectively. The matrix is assumed
1862 to have n columns and an arbitrary number of rows.
1863 On entry, len(*) must contain the length of each column.
1864 On exit, a(*) and inum(*) are sorted, jnum(*) = 0, and
1865 loc(j) points to the start of column j.
1866 lu1or2 is derived from mc20ad, a routine in the Harwell
1867 Subroutine Library, author J. K. Reid.
1868 ------------------------------------------------------------------
1869 xx Feb 1985: Original version.
1870 17 Oct 2000: a, inum, jnum now have size lena to allow nelem = 0.
1871 ================================================================== */
1872 void LU1OR2(LUSOLrec *LUSOL)
1873 {
1874 REAL ACE, ACEP;
1875 int L, J, I, JCE, ICE, ICEP, JCEP, JA, JB;
1876
1877 /* Set loc(j) to point to the beginning of column j. */
1878 L = 1;
1879 for(J = 1; J <= LUSOL->n; J++) {
1880 LUSOL->locc[J] = L;
1881 L += LUSOL->lenc[J];
1882 }
1883 /* Sort the elements into column order.
1884 The algorithm is an in-place sort and is of order numa. */
1885 for(I = 1; I <= LUSOL->nelem; I++) {
1886 /* Establish the current entry. */
1887 JCE = LUSOL->indr[I];
1888 if(JCE==0)
1889 continue;
1890 ACE = LUSOL->a[I];
1891 ICE = LUSOL->indc[I];
1892 LUSOL->indr[I] = 0;
1893 /* Chain from current entry. */
1894 for(J = 1; J <= LUSOL->nelem; J++) {
1895 /* The current entry is not in the correct position.
1896 Determine where to store it. */
1897 L = LUSOL->locc[JCE];
1898 LUSOL->locc[JCE]++;
1899 /* Save the contents of that location. */
1900 ACEP = LUSOL->a[L];
1901 ICEP = LUSOL->indc[L];
1902 JCEP = LUSOL->indr[L];
1903 /* Store current entry. */
1904 LUSOL->a[L] = ACE;
1905 LUSOL->indc[L] = ICE;
1906 LUSOL->indr[L] = 0;
1907 /* If next current entry needs to be processed,
1908 copy it into current entry. */
1909 if(JCEP==0)
1910 break;
1911 ACE = ACEP;
1912 ICE = ICEP;
1913 JCE = JCEP;
1914 }
1915 }
1916 /* Reset loc(j) to point to the start of column j. */
1917 JA = 1;
1918 for(J = 1; J <= LUSOL->n; J++) {
1919 JB = LUSOL->locc[J];
1920 LUSOL->locc[J] = JA;
1921 JA = JB;
1922 }
1923 }
1924
1925 /* ==================================================================
1926 lu1or3 looks for duplicate elements in an m by n matrix A
1927 defined by the column list indc, lenc, locc.
1928 iw is used as a work vector of length m.
1929 ------------------------------------------------------------------
1930 xx Feb 1985: Original version.
1931 17 Oct 2000: indc, indr now have size lena to allow nelem = 0.
1932 ================================================================== */
1933 void LU1OR3(LUSOLrec *LUSOL, int *LERR, int *INFORM)
1934 {
1935 int I, J, L1, L2, L;
1936
1937 #ifdef LUSOLFastClear
1938 MEMCLEAR((LUSOL->ip+1), LUSOL->m);
1939 #else
1940 for(I = 1; I <= LUSOL->m; I++)
1941 LUSOL->ip[I] = ZERO;
1942 #endif
1943
1944 for(J = 1; J <= LUSOL->n; J++) {
1945 if(LUSOL->lenc[J]>0) {
1946 L1 = LUSOL->locc[J];
1947 L2 = (L1+LUSOL->lenc[J])-1;
1948 for(L = L1; L <= L2; L++) {
1949 I = LUSOL->indc[L];
1950 if(LUSOL->ip[I]==J)
1951 goto x910;
1952 LUSOL->ip[I] = J;
1953 }
1954 }
1955 }
1956 *INFORM = LUSOL_INFORM_LUSUCCESS;
1957 return;
1958 x910:
1959 *LERR = L;
1960 *INFORM = LUSOL_INFORM_LUSINGULAR;
1961 }
1962
1963 /* ==================================================================
1964 lu1or4 constructs a row list indr, locr
1965 from a corresponding column list indc, locc,
1966 given the lengths of both columns and rows in lenc, lenr.
1967 ------------------------------------------------------------------
1968 xx Feb 1985: Original version.
1969 17 Oct 2000: indc, indr now have size lena to allow nelem = 0.
1970 ================================================================== */
1971 void LU1OR4(LUSOLrec *LUSOL)
1972 {
1973 int L, I, L2, J, JDUMMY, L1, LR;
1974
1975 /* Initialize locr(i) to point just beyond where the
1976 last component of row i will be stored. */
1977 L = 1;
1978 for(I = 1; I <= LUSOL->m; I++) {
1979 L += LUSOL->lenr[I];
1980 LUSOL->locr[I] = L;
1981 }
1982 /* By processing the columns backwards and decreasing locr(i)
1983 each time it is accessed, it will end up pointing to the
1984 beginning of row i as required. */
1985 L2 = LUSOL->nelem;
1986 J = LUSOL->n+1;
1987 for(JDUMMY = 1; JDUMMY <= LUSOL->n; JDUMMY++) {
1988 J = J-1;
1989 if(LUSOL->lenc[J]>0) {
1990 L1 = LUSOL->locc[J];
1991 for(L = L1; L <= L2; L++) {
1992 I = LUSOL->indc[L];
1993 LR = LUSOL->locr[I]-1;
1994 LUSOL->locr[I] = LR;
1995 LUSOL->indr[LR] = J;
1996 }
1997 L2 = L1-1;
1998 }
1999 }
2000 }
2001
2002 /* ==================================================================
2003 lu1pen deals with pending fill-in in the row file.
2004 ------------------------------------------------------------------
2005 ifill(ll) says if a row involved in the new column of L
2006 has to be updated. If positive, it is the total
2007 length of the final updated row.
2008 jfill(lu) says if a column involved in the new row of U
2009 contains any pending fill-ins. If positive, it points
2010 to the first fill-in in the column that has yet to be
2011 added to the row file.
2012 ------------------------------------------------------------------
2013 16 Apr 1989: First version of lu1pen.
2014 23 Mar 2001: ilast used and updated.
2015 ================================================================== */
2016 void LU1PEN(LUSOLrec *LUSOL, int NSPARE, int *ILAST,
2017 int LPIVC1, int LPIVC2, int LPIVR1, int LPIVR2,
2018 int *LROW, int IFILL[], int JFILL[])
2019 {
2020 int LL, LC, L, I, LR1, LR2, LR, LU, J, LC1, LC2, LAST;
2021
2022 LL = 0;
2023 for(LC = LPIVC1; LC <= LPIVC2; LC++) {
2024 LL++;
2025 if(IFILL[LL]==0)
2026 continue;
2027 /* Another row has pending fill.
2028 First, add some spare space at the }
2029 of the current last row. */
2030 #if 1
2031 LC1 = (*LROW)+1;
2032 LC2 = (*LROW)+NSPARE;
2033 *LROW = LC2;
2034 for(L = LC1; L <= LC2; L++) {
2035 #else
2036 for(L = (*LROW)+1; L <= (*LROW)+NSPARE; L++) {
2037 *LROW = L; /* ******* ERROR ???? */
2038 #endif
2039 LUSOL->indr[L] = 0;
2040 }
2041 /* Now move row i to the end of the row file. */
2042 I = LUSOL->indc[LC];
2043 *ILAST = I;
2044 LR1 = LUSOL->locr[I];
2045 LR2 = (LR1+LUSOL->lenr[I])-1;
2046 LUSOL->locr[I] = (*LROW)+1;
2047 for(LR = LR1; LR <= LR2; LR++) {
2048 (*LROW)++;
2049 LUSOL->indr[*LROW] = LUSOL->indr[LR];
2050 LUSOL->indr[LR] = 0;
2051 }
2052 (*LROW) += IFILL[LL];
2053 }
2054 /* Scan all columns of D and insert the pending fill-in
2055 into the row file. */
2056 LU = 1;
2057 for(LR = LPIVR1; LR <= LPIVR2; LR++) {
2058 LU++;
2059 if(JFILL[LU]==0)
2060 continue;
2061 J = LUSOL->indr[LR];
2062 LC1 = (LUSOL->locc[J]+JFILL[LU])-1;
2063 LC2 = (LUSOL->locc[J]+LUSOL->lenc[J])-1;
2064 for(LC = LC1; LC <= LC2; LC++) {
2065 I = LUSOL->indc[LC]-LUSOL->m;
2066 if(I>0) {
2067 LUSOL->indc[LC] = I;
2068 LAST = LUSOL->locr[I]+LUSOL->lenr[I];
2069 LUSOL->indr[LAST] = J;
2070 LUSOL->lenr[I]++;
2071 }
2072 }
2073 }
2074 }
2075
2076
2077 /* ==================================================================
2078 lu1fad is a driver for the numerical phase of lu1fac.
2079 At each stage it computes a column of L and a row of U,
2080 using a Markowitz criterion to select the pivot element,
2081 subject to a stability criterion that bounds the elements of L.
2082 ------------------------------------------------------------------
2083 Local variables
2084 ---------------
2085 lcol is the length of the column file. It points to the last
2086 nonzero in the column list.
2087 lrow is the analogous quantity for the row file.
2088 lfile is the file length (lcol or lrow) after the most recent
2089 compression of the column list or row list.
2090 nrowd and ncold are the number of rows and columns in the
2091 matrix defined by the pivot column and row. They are the
2092 dimensions of the submatrix D being altered at this stage.
2093 melim and nelim are the number of rows and columns in the
2094 same matrix D, excluding the pivot column and row.
2095 mleft and nleft are the number of rows and columns
2096 still left to be factored.
2097 nzchng is the increase in nonzeros in the matrix that remains
2098 to be factored after the current elimination
2099 (usually negative).
2100 nzleft is the number of nonzeros still left to be factored.
2101 nspare is the space we leave at the end of the last row or
2102 column whenever a row or column is being moved to the }
2103 of its file. nspare = 1 or 2 might help reduce the
2104 number of file compressions when storage is tight.
2105 The row and column ordering permutes A into the form
2106 ------------------------
2107 \ |
2108 \ U1 |
2109 \ |
2110 --------------------
2111 |\
2112 | \
2113 | \
2114 P A Q = | \
2115 | \
2116 | --------------
2117 | | |
2118 | | |
2119 | L1 | A2 |
2120 | | |
2121 | | |
2122 --------------------
2123 where the block A2 is factored as A2 = L2 U2.
2124 The phases of the factorization are as follows.
2125 Utri is true when U1 is being determined.
2126 Any column of length 1 is accepted immediately (if TPP).
2127 Ltri is true when L1 is being determined.
2128 lu1mar exits as soon as an acceptable pivot is found
2129 in a row of length 1.
2130 spars1 is true while the density of the (modified) A2 is less
2131 than the parameter dens1 = parmlu(7) = 0.3 say.
2132 lu1mar searches maxcol columns and maxrow rows,
2133 where maxcol = luparm(3), maxrow = maxcol - 1.
2134 lu1mxc is used to keep the biggest element at the top
2135 of all remaining columns.
2136 spars2 is true while the density of the modified A2 is less
2137 than the parameter dens2 = parmlu(8) = 0.6 say.
2138 lu1mar searches maxcol columns and no rows.
2139 lu1mxc could fix up only the first maxcol cols (with TPP).
2140 22 Sep 2000: For simplicity, lu1mxc fixes all
2141 modified cols.
2142 dense is true once the density of A2 reaches dens2.
2143 lu1mar searches only 1 column (the shortest).
2144 lu1mxc could fix up only the first column (with TPP).
2145 ------------------------------------------------------------------
2146 00 Jan 1986 Version documented in LUSOL paper:
2147 Gill, Murray, Saunders and Wright (1987),
2148 Maintaining LU factors of a general sparse matrix,
2149 Linear algebra and its applications 88/89, 239-270.
2150 02 Feb 1989 Following Suhl and Aittoniemi (1987), the largest
2151 element in each column is now kept at the start of
2152 the column, i.e. in position locc(j) of a and indc.
2153 This should speed up the Markowitz searches.
2154 To save time on highly triangular matrices, we wait
2155 until there are no further columns of length 1
2156 before setting and maintaining that property.
2157 12 Apr 1989 ipinv and iqinv added (inverses of ip and iq)
2158 to save searching ip and iq for rows and columns
2159 altered in each elimination step. (Used in lu1pq2)
2160 19 Apr 1989 Code segmented to reduce its size.
2161 lu1gau does most of the Gaussian elimination work.
2162 lu1mar does just the Markowitz search.
2163 lu1mxc moves biggest elements to top of columns.
2164 lu1pen deals with pending fill-in in the row list.
2165 lu1pq2 updates the row and column permutations.
2166 26 Apr 1989 maxtie replaced by maxcol, maxrow in the Markowitz
2167 search. maxcol, maxrow change as density increases.
2168 25 Oct 1993 keepLU implemented.
2169 07 Feb 1994 Exit main loop early to finish off with a dense LU.
2170 densLU tells lu1fad whether to do it.
2171 21 Dec 1994 Bug fixed. nrank was wrong after the call to lu1ful.
2172 12 Nov 1999 A parallel version of dcopy gave trouble in lu1ful
2173 during left-shift of dense matrix D within a(*).
2174 Fixed this unexpected problem here in lu1fad
2175 by making sure the first and second D don't overlap.
2176 13 Sep 2000 TCP (Threshold Complete Pivoting) implemented.
2177 lu2max added
2178 (finds aijmax from biggest elems in each col).
2179 Utri, Ltri and Spars1 phases apply.
2180 No switch to Dense CP yet. (Only TPP switches.)
2181 14 Sep 2000 imax needed to remember row containing aijmax.
2182 22 Sep 2000 For simplicity, lu1mxc always fixes all modified cols.
2183 (TPP spars2 used to fix just the first maxcol cols.)
2184 08 Nov 2000: Speed up search for aijmax.
2185 Don't need to search all columns if the elimination
2186 didn't alter the col containing the current aijmax.
2187 21 Nov 2000: lu1slk implemented for Utri phase with TCP
2188 to guard against deceptive triangular matrices.
2189 (Utri used to have aijtol >= 0.9999 to include
2190 slacks, but this allows other 1s to be accepted.)
2191 Utri now accepts slacks, but applies normal aijtol
2192 test to other pivots.
2193 28 Nov 2000: TCP with empty cols must call lu1mxc and lu2max
2194 with ( lq1, n, ... ), not just ( 1, n, ... ).
2195 23 Mar 2001: lu1fad bug with TCP.
2196 A col of length 1 might not be accepted as a pivot.
2197 Later it appears in a pivot row and temporarily
2198 has length 0 (when pivot row is removed
2199 but before the column is filled in). If it is the
2200 last column in storage, the preceding col also thinks
2201 it is "last". Trouble arises when the preceding col
2202 needs fill-in -- it overlaps the real "last" column.
2203 (Very rarely, same trouble might have happened if
2204 the drop tolerance caused columns to have length 0.)
2205 Introduced ilast to record the last row in row file,
2206 jlast to record the last col in col file.
2207 lu1rec returns ilast = indr(lrow + 1)
2208 or jlast = indc(lcol + 1).
2209 (Should be an output parameter, but didn't want to
2210 alter lu1rec's parameter list.)
2211 lu1rec also treats empty rows or cols safely.
2212 (Doesn't eliminate them!)
2213 26 Apr 2002: Heap routines added for TCP.
2214 lu2max no longer needed.
2215 imax, jmax used only for printing.
2216 01 May 2002: lu1DCP implemented (dense complete pivoting).
2217 Both TPP and TCP now switch to dense LU
2218 when density exceeds dens2.
2219 06 May 2002: In dense mode, store diag(U) in natural order.
2220 09 May 2002: lu1mCP implemented (Markowitz TCP via heap).
2221 11 Jun 2002: lu1mRP implemented (Markowitz TRP).
2222 28 Jun 2002: Fixed call to lu1mxr.
2223 14 Dec 2002: lu1mSP implemented (Markowitz TSP).
2224 15 Dec 2002: Both TPP and TSP can grab cols of length 1
2225 during Utri.
2226 ================================================================== */
2227 void LU1FAD(LUSOLrec *LUSOL,
2228 #ifdef ClassicHamaxR
2229 int LENA2, int LENH, REAL HA[], int HJ[], int HK[], REAL AMAXR[],
2230 #endif
2231 int *INFORM, int *LENL, int *LENU, int *MINLEN,
2232 int *MERSUM, int *NUTRI, int *NLTRI,
2233 int *NDENS1, int *NDENS2, int *NRANK,
2234 REAL *LMAX, REAL *UMAX, REAL *DUMAX, REAL *DUMIN, REAL *AKMAX)
2235 {
2236 MYBOOL UTRI, LTRI, SPARS1, SPARS2, DENSE, DENSLU, KEEPLU, TCP, TPP, TRP,TSP;
2237 int HLEN, HOPS, H, LPIV, LPRINT, MAXCOL, MAXROW, ILAST, JLAST, LFILE, LROW, LCOL,
2238 MINMN, MAXMN, NZLEFT, NSPARE, LU1, KK, J, LC, MLEFT, NLEFT, NROWU,
2239 LQ1, LQ2, JBEST, LQ, I, IBEST, MBEST, LEND, NFREE, LD, NCOLD, NROWD,
2240 MELIM, NELIM, JMAX, IMAX, LL1, LSAVE, LFREE, LIMIT, MINFRE, LPIVR, LPIVR1, LPIVR2,
2241 L, LPIVC, LPIVC1, LPIVC2, KBEST, LU, LR, LENJ, LC1, LAST, LL, LS,
2242 LENI, LR1, LFIRST, NFILL, NZCHNG, K, MRANK, NSING;
2243 REAL LIJ, LTOL, SMALL, USPACE, DENS1, DENS2, AIJMAX, AIJTOL, AMAX, ABEST, DIAG, V;
2244 #ifdef ClassicHamaxR
2245 int LDIAGU;
2246 #else
2247 int LENA2 = LUSOL->lena;
2248 #endif
2249
2250 #ifdef UseTimer
2251 int eltime, mktime, ntime;
2252 timer ( "start", 3 );
2253 ntime = LUSOL->n / 4;
2254 #endif
2255
2256 #ifdef ForceInitialization
2257 AIJMAX = 0;
2258 AIJTOL = 0;
2259 HLEN = 0;
2260 JBEST = 0;
2261 IBEST = 0;
2262 MBEST = 0;
2263 LEND = 0;
2264 LD = 0;
2265 #endif
2266
2267 LPRINT = LUSOL->luparm[LUSOL_IP_PRINTLEVEL];
2268 MAXCOL = LUSOL->luparm[LUSOL_IP_MARKOWITZ_MAXCOL];
2269 LPIV = LUSOL->luparm[LUSOL_IP_PIVOTTYPE];
2270 KEEPLU = (MYBOOL) (LUSOL->luparm[LUSOL_IP_KEEPLU]!=FALSE);
2271 /* Threshold Partial Pivoting (normal). */
2272 TPP = (MYBOOL) (LPIV==LUSOL_PIVMOD_TPP);
2273 /* Threshold Rook Pivoting */
2274 TRP = (MYBOOL) (LPIV==LUSOL_PIVMOD_TRP);
2275 /* Threshold Complete Pivoting. */
2276 TCP = (MYBOOL) (LPIV==LUSOL_PIVMOD_TCP);
2277 /* Threshold Symmetric Pivoting. */
2278 TSP = (MYBOOL) (LPIV==LUSOL_PIVMOD_TSP);
2279 DENSLU = FALSE;
2280 MAXROW = MAXCOL-1;
2281 /* Assume row m is last in the row file. */
2282 ILAST = LUSOL->m;
2283 /* Assume col n is last in the col file. */
2284 JLAST = LUSOL->n;
2285 LFILE = LUSOL->nelem;
2286 LROW = LUSOL->nelem;
2287 LCOL = LUSOL->nelem;
2288 MINMN = MIN(LUSOL->m,LUSOL->n);
2289 MAXMN = MAX(LUSOL->m,LUSOL->n);
2290 NZLEFT = LUSOL->nelem;
2291 NSPARE = 1;
2292
2293 if(KEEPLU)
2294 LU1 = LENA2+1;
2295 else {
2296 /* Store only the diagonals of U in the top of memory. */
2297 #ifdef ClassicdiagU
2298 LDIAGU = LENA2-LUSOL->n;
2299 LU1 = LDIAGU+1;
2300 LUSOL->diagU = LUSOL->a+LDIAGU;
2301 #else
2302 LU1 = LENA2+1;
2303 #endif
2304 }
2305
2306 LTOL = LUSOL->parmlu[LUSOL_RP_FACTORMAX_Lij];
2307 SMALL = LUSOL->parmlu[LUSOL_RP_ZEROTOLERANCE];
2308 USPACE = LUSOL->parmlu[LUSOL_RP_COMPSPACE_U];
2309 DENS1 = LUSOL->parmlu[LUSOL_RP_MARKOWITZ_CONLY];
2310 DENS2 = LUSOL->parmlu[LUSOL_RP_MARKOWITZ_DENSE];
2311 UTRI = TRUE;
2312 LTRI = FALSE;
2313 SPARS1 = FALSE;
2314 SPARS2 = FALSE;
2315 DENSE = FALSE;
2316 /* Check parameters. */
2317 SETMAX(LTOL,1.0001E+0);
2318 SETMIN(DENS1,DENS2);
2319 /* Initialize output parameters.
2320 lenL, lenU, minlen, mersum, nUtri, nLtri, ndens1, ndens2, nrank
2321 are already initialized by lu1fac. */
2322 *LMAX = ZERO;
2323 *UMAX = ZERO;
2324 *DUMAX = ZERO;
2325 *DUMIN = LUSOL_BIGNUM;
2326 if(LUSOL->nelem==0)
2327 *DUMIN = ZERO;
2328 *AKMAX = ZERO;
2329 HOPS = 0;
2330 /* More initialization.
2331 Don't worry yet about lu1mxc. */
2332 if(TPP || TSP) {
2333 AIJMAX = ZERO;
2334 AIJTOL = ZERO;
2335 HLEN = 1;
2336 /* TRP or TCP */
2337 }
2338 else {
2339 /* Move biggest element to top of each column.
2340 Set w(*) to mark slack columns (unit vectors). */
2341 LU1MXC(LUSOL, 1,LUSOL->n,LUSOL->iq);
2342 LU1SLK(LUSOL);
2343 }
2344 if(TRP)
2345 /* Find biggest element in each row. */
2346 #ifdef ClassicHamaxR
2347 LU1MXR(LUSOL, 1,LUSOL->m,LUSOL->ip,AMAXR);
2348 #else
2349 LU1MXR(LUSOL, 1,LUSOL->m,LUSOL->ip,LUSOL->amaxr);
2350 #endif
2351
2352 if(TCP) {
2353 /* Set Ha(1:Hlen) = biggest element in each column,
2354 Hj(1:Hlen) = corresponding column indices. */
2355 HLEN = 0;
2356 for(KK = 1; KK <= LUSOL->n; KK++) {
2357 HLEN++;
2358 J = LUSOL->iq[KK];
2359 LC = LUSOL->locc[J];
2360 #ifdef ClassicHamaxR
2361 HA[HLEN] = fabs(LUSOL->a[LC]);
2362 HJ[HLEN] = J;
2363 HK[J] = HLEN;
2364 #else
2365 LUSOL->Ha[HLEN] = fabs(LUSOL->a[LC]);
2366 LUSOL->Hj[HLEN] = J;
2367 LUSOL->Hk[J] = HLEN;
2368 #endif
2369 }
2370 /* Build the heap, creating new Ha, Hj and setting Hk(1:Hlen). */
2371 #ifdef ClassicHamaxR
2372 HBUILD(HA,HJ,HK,HLEN,&HOPS);
2373 #else
2374 HBUILD(LUSOL->Ha,LUSOL->Hj,LUSOL->Hk,HLEN,&HOPS);
2375 #endif
2376 }
2377 /* ------------------------------------------------------------------
2378 Start of main loop.
2379 ------------------------------------------------------------------ */
2380 MLEFT = LUSOL->m+1;
2381 NLEFT = LUSOL->n+1;
2382 for(NROWU = 1; NROWU <= MINMN; NROWU++) {
2383 #ifdef UseTimer
2384 mktime = (nrowu / ntime) + 4;
2385 eltime = (nrowu / ntime) + 9;
2386 #endif
2387 MLEFT--;
2388 NLEFT--;
2389 /* Bail out if there are no nonzero rows left. */
2390 if(LUSOL->iploc[1]>LUSOL->m)
2391 goto x900;
2392 /* For TCP, the largest Aij is at the top of the heap. */
2393 if(TCP) {
2394 /*
2395 Marvelously easy */
2396 #ifdef ClassicHamaxR
2397 AIJMAX = HA[1];
2398 #else
2399 AIJMAX = LUSOL->Ha[1];
2400 #endif
2401 SETMAX(*AKMAX,AIJMAX);
2402 AIJTOL = AIJMAX/LTOL;
2403 }
2404 /* ===============================================================
2405 Find a suitable pivot element.
2406 =============================================================== */
2407 if(UTRI) {
2408 /* ------------------------------------------------------------
2409 So far all columns have had length 1.
2410 We are still looking for the (backward) triangular part of A
2411 that forms the first rows and columns of U.
2412 ------------------------------------------------------------ */
2413 LQ1 = LUSOL->iqloc[1];
2414 LQ2 = LUSOL->n;
2415 if(LUSOL->m>1)
2416 LQ2 = LUSOL->iqloc[2]-1;
2417 /* There are more cols of length 1. */
2418 if(LQ1<=LQ2) {
2419 if(TPP || TSP) {
2420 /* Grab the first one. */
2421 JBEST = LUSOL->iq[LQ1];
2422 /* Scan all columns of length 1 ... TRP or TCP */
2423 }
2424 else {
2425 JBEST = 0;
2426 for(LQ = LQ1; LQ <= LQ2; LQ++) {
2427 J = LUSOL->iq[LQ];
2428 /* Accept a slack */
2429 if(LUSOL->w[J]>ZERO) {
2430 JBEST = J;
2431 goto x250;
2432 }
2433 LC = LUSOL->locc[J];
2434 AMAX = fabs(LUSOL->a[LC]);
2435 if(TRP) {
2436 I = LUSOL->indc[LC];
2437 #ifdef ClassicHamaxR
2438 AIJTOL = AMAXR[I]/LTOL;
2439 #else
2440 AIJTOL = LUSOL->amaxr[I]/LTOL;
2441 #endif
2442 }
2443 if(AMAX>=AIJTOL) {
2444 JBEST = J;
2445 goto x250;
2446 }
2447 }
2448 }
2449 x250:
2450 if(JBEST>0) {
2451 LC = LUSOL->locc[JBEST];
2452 IBEST = LUSOL->indc[LC];
2453 MBEST = 0;
2454 goto x300;
2455 }
2456 }
2457 /* This is the end of the U triangle.
2458 We will not return to this part of the code.
2459 TPP and TSP call lu1mxc for the first time
2460 (to move biggest element to top of each column). */
2461 if(LPRINT>=LUSOL_MSG_PIVOT)
2462 LUSOL_report(LUSOL, 0, "Utri ended. spars1 = TRUE\n");
2463 UTRI = FALSE;
2464 LTRI = TRUE;
2465 SPARS1 = TRUE;
2466 *NUTRI = NROWU-1;
2467 if(TPP || TSP)
2468 LU1MXC(LUSOL, LQ1,LUSOL->n,LUSOL->iq);
2469 }
2470 if(SPARS1) {
2471 /* ------------------------------------------------------------
2472 Perform a Markowitz search.
2473 Search cols of length 1, then rows of length 1,
2474 then cols of length 2, then rows of length 2, etc.
2475 ------------------------------------------------------------ */
2476 #ifdef UseTimer
2477 timer ( "start", mktime );
2478 #endif
2479 /* 12 Jun 2002: Next line disables lu1mCP below
2480 if (TPP) then */
2481 if(TPP || TCP) {
2482 LU1MAR(LUSOL, MAXMN,TCP,AIJTOL,LTOL,MAXCOL,MAXROW,&IBEST,&JBEST,&MBEST);
2483 }
2484 else if(TRP) {
2485 #ifdef ClassicHamaxR
2486 LU1MRP(LUSOL, MAXMN,LTOL,MAXCOL,MAXROW,&IBEST,&JBEST,&MBEST,AMAXR);
2487 #else
2488 LU1MRP(LUSOL, MAXMN,LTOL,MAXCOL,MAXROW,&IBEST,&JBEST,&MBEST,LUSOL->amaxr);
2489 #endif
2490 /* else if (TCP) {
2491 lu1mCP( m , n , lena , aijtol,
2492 ibest, jbest , mbest ,
2493 a , indc , indr ,
2494 lenc , lenr , locc ,
2495 Hlen , Ha , Hj ) */
2496 }
2497 else if(TSP) {
2498 LU1MSP(LUSOL, MAXMN,LTOL,MAXCOL,&IBEST,&JBEST,&MBEST);
2499 if(IBEST==0)
2500 goto x990;
2501 }
2502 #ifdef UseTimer
2503 timer ( "finish", mktime );
2504 #endif
2505 if(LTRI) {
2506 /* So far all rows have had length 1.
2507 We are still looking for the (forward) triangle of A
2508 that forms the first rows and columns of L. */
2509 if(MBEST>0) {
2510 LTRI = FALSE;
2511 *NLTRI = NROWU-1-*NUTRI;
2512 if(LPRINT>=LUSOL_MSG_PIVOT)
2513 LUSOL_report(LUSOL, 0, "Ltri ended.\n");
2514 }
2515 }
2516 else {
2517 /* See if what's left is as dense as dens1. */
2518 if(NZLEFT>=(DENS1*MLEFT)*NLEFT) {
2519 SPARS1 = FALSE;
2520 SPARS2 = TRUE;
2521 *NDENS1 = NLEFT;
2522 MAXROW = 0;
2523 if(LPRINT>=LUSOL_MSG_PIVOT)
2524 LUSOL_report(LUSOL, 0, "spars1 ended. spars2 = TRUE\n");
2525 }
2526 }
2527 }
2528 else if(SPARS2 || DENSE) {
2529 /* ------------------------------------------------------------
2530 Perform a restricted Markowitz search,
2531 looking at only the first maxcol columns. (maxrow = 0.)
2532 ------------------------------------------------------------ */
2533 #ifdef UseTimer
2534 timer ( "start", mktime );
2535 #endif
2536 /* 12 Jun 2002: Next line disables lu1mCP below
2537 if (TPP) then */
2538 if(TPP || TCP) {
2539 LU1MAR(LUSOL, MAXMN,TCP,AIJTOL,LTOL,MAXCOL,MAXROW,&IBEST,&JBEST,&MBEST);
2540 }
2541 else if(TRP) {
2542 #ifdef ClassicHamaxR
2543 LU1MRP(LUSOL, MAXMN,LTOL,MAXCOL,MAXROW,&IBEST,&JBEST,&MBEST,AMAXR);
2544 #else
2545 LU1MRP(LUSOL, MAXMN,LTOL,MAXCOL,MAXROW,&IBEST,&JBEST,&MBEST,LUSOL->amaxr);
2546 #endif
2547 /* else if (TCP) {
2548 lu1mCP( m , n , lena , aijtol,
2549 ibest, jbest , mbest ,
2550 a , indc , indr ,
2551 lenc , lenr , locc ,
2552 Hlen , Ha , Hj ) */
2553 }
2554 else if(TSP) {
2555 LU1MSP(LUSOL, MAXMN,LTOL,MAXCOL,&IBEST,&JBEST,&MBEST);
2556 if(IBEST==0)
2557 goto x985;
2558 }
2559 #ifdef UseTimer
2560 timer ( "finish", mktime );
2561 #endif
2562 /* See if what's left is as dense as dens2. */
2563 if(SPARS2) {
2564 if(NZLEFT>=(DENS2*MLEFT)*NLEFT) {
2565 SPARS2 = FALSE;
2566 DENSE = TRUE;
2567 *NDENS2 = NLEFT;
2568 MAXCOL = 1;
2569 if(LPRINT>=LUSOL_MSG_PIVOT)
2570 LUSOL_report(LUSOL, 0, "spars2 ended. dense = TRUE\n");
2571 }
2572 }
2573 }
2574 /* ---------------------------------------------------------------
2575 See if we can finish quickly.
2576 --------------------------------------------------------------- */
2577 if(DENSE) {
2578 LEND = MLEFT*NLEFT;
2579 NFREE = LU1-1;
2580 if(NFREE>=2*LEND) {
2581 /* There is room to treat the remaining matrix as
2582 a dense matrix D.
2583 We may have to compress the column file first.
2584 12 Nov 1999: D used to be put at the
2585 beginning of free storage (lD = lcol + 1).
2586 Now put it at the end (lD = lu1 - lenD)
2587 so the left-shift in lu1ful will not
2588 involve overlapping storage
2589 (fatal with parallel dcopy).
2590 */
2591 DENSLU = TRUE;
2592 *NDENS2 = NLEFT;
2593 LD = LU1-LEND;
2594 if(LCOL>=LD) {
2595 LU1REC(LUSOL, LUSOL->n,TRUE,&LCOL,
2596 LUSOL->indc,LUSOL->lenc,LUSOL->locc);
2597 LFILE = LCOL;
2598 JLAST = LUSOL->indc[LCOL+1];
2599 }
2600 goto x900;
2601 }
2602 }
2603 /* ===============================================================
2604 The best aij has been found.
2605 The pivot row ibest and the pivot column jbest
2606 Define a dense matrix D of size nrowd by ncold.
2607 =============================================================== */
2608 x300:
2609 NCOLD = LUSOL->lenr[IBEST];
2610 NROWD = LUSOL->lenc[JBEST];
2611 MELIM = NROWD-1;
2612 NELIM = NCOLD-1;
2613 (*MERSUM) += MBEST;
2614 (*LENL) += MELIM;
2615 (*LENU) += NCOLD;
2616 if(LPRINT>=LUSOL_MSG_PIVOT) {
2617 if(NROWU==1)
2618 LUSOL_report(LUSOL, 0, "lu1fad debug:\n");
2619 if(TPP || TRP || TSP) {
2620 LUSOL_report(LUSOL, 0, "nrowu:%7d i,jbest:%7d,%7d nrowd,ncold:%6d,%6d\n",
2621 NROWU, IBEST,JBEST, NROWD,NCOLD);
2622 /* TCP */
2623 }
2624 else {
2625 #ifdef ClassicHamaxR
2626 JMAX = HJ[1];
2627 #else
2628 JMAX = LUSOL->Hj[1];
2629 #endif
2630 IMAX = LUSOL->indc[LUSOL->locc[JMAX]];
2631 LUSOL_report(LUSOL, 0, "nrowu:%7d i,jbest:%7d,%7d nrowd,ncold:%6d,%6d i,jmax:%7d,%7d aijmax:%g\n",
2632 NROWU, IBEST,JBEST, NROWD,NCOLD, IMAX,JMAX, AIJMAX);
2633 }
2634 }
2635 /* ===============================================================
2636 Allocate storage for the next column of L and next row of U.
2637 Initially the top of a, indc, indr are used as follows:
2638 ncold melim ncold melim
2639 a |...........|...........|ujbest..ujn|li1......lim|
2640 indc |...........| lenr(i) | lenc(j) | markl(i) |
2641 indr |...........| iqloc(i) | jfill(j) | ifill(i) |
2642 ^ ^ ^ ^ ^
2643 lfree lsave lu1 ll1 oldlu1
2644 Later the correct indices are inserted:
2645 indc | | | |i1........im|
2646 indr | | |jbest....jn|ibest..ibest|
2647 =============================================================== */
2648 if(!KEEPLU) {
2649 /* Always point to the top spot.
2650 Only the current column of L and row of U will
2651 take up space, overwriting the previous ones. */
2652 #ifdef ClassicHamaxR
2653 LU1 = LDIAGU+1;
2654 #else
2655 LU1 = LENA2+1;
2656 #endif
2657 }
2658 /* Update (left-shift) pointers to make room for the new data */
2659 LL1 = LU1-MELIM;
2660 LU1 = LL1-NCOLD;
2661 LSAVE = LU1-NROWD;
2662 LFREE = LSAVE-NCOLD;
2663
2664 /* Check if we need to allocate more memory, and allocate if necessary */
2665 #if 0 /* Proposal by Michael A. Saunders (logic based on Markowitz' rule) */
2666 L = NROWD*NCOLD;
2667
2668 /* Try to avoid future expansions by anticipating further updates - KE extension */
2669 if(LUSOL->luparm[LUSOL_IP_UPDATELIMIT] > 0)
2670 #if 1
2671 L *= (int) (log(LUSOL->luparm[LUSOL_IP_UPDATELIMIT]-LUSOL->luparm[LUSOL_IP_UPDATECOUNT]+2.0) + 1);
2672 #else
2673 L *= (LUSOL->luparm[LUSOL_IP_UPDATELIMIT]-LUSOL->luparm[LUSOL_IP_UPDATECOUNT]) / 2 + 1;
2674 #endif
2675
2676 #else /* Version by Kjell Eikland (from luparm[LUSOL_IP_MINIMUMLENA] and safety margin) */
2677 L = (KEEPLU ? MAX(LROW, LCOL) + 2*(LUSOL->m+LUSOL->n) : 0);
2678 L *= LUSOL_MULT_nz_a;
2679 SETMAX(L, NROWD*NCOLD);
2680 #endif
2681
2682 /* Do the memory expansion */
2683 if((L > LFREE-LCOL) && LUSOL_expand_a(LUSOL, &L, &LFREE)) {
2684 LL1 += L;
2685 LU1 += L;
2686 LSAVE += L;
2687 #ifdef ClassicdiagU
2688 LUSOL->diagU += L;
2689 #endif
2690 #ifdef ClassicHamaxR
2691 HA += L;
2692 HJ += L;
2693 HK += L;
2694 AMAXR += L;
2695 #endif
2696 }
2697 LIMIT = (int) (USPACE*LFILE)+LUSOL->m+LUSOL->n+1000;
2698
2699 /* Make sure the column file has room.
2700 Also force a compression if its length exceeds a certain limit. */
2701 #ifdef StaticMemAlloc
2702 MINFRE = NCOLD+MELIM;
2703 #else
2704 MINFRE = NROWD*NCOLD;
2705 #endif
2706 NFREE = LFREE-LCOL;
2707 if(NFREE<MINFRE || LCOL>LIMIT) {
2708 LU1REC(LUSOL, LUSOL->n,TRUE,&LCOL,
2709 LUSOL->indc,LUSOL->lenc,LUSOL->locc);
2710 LFILE = LCOL;
2711 JLAST = LUSOL->indc[LCOL+1];
2712 NFREE = LFREE-LCOL;
2713 if(NFREE<MINFRE)
2714 goto x970;
2715 }
2716 /* Make sure the row file has room. */
2717 #ifdef StaticMemAlloc
2718 MINFRE = NCOLD+MELIM;
2719 #else
2720 MINFRE = NROWD*NCOLD;
2721 #endif
2722 NFREE = LFREE-LROW;
2723 if(NFREE<MINFRE || LROW>LIMIT) {
2724 LU1REC(LUSOL, LUSOL->m,FALSE,&LROW,
2725 LUSOL->indr,LUSOL->lenr,LUSOL->locr);
2726 LFILE = LROW;
2727 ILAST = LUSOL->indr[LROW+1];
2728 NFREE = LFREE-LROW;
2729 if(NFREE<MINFRE)
2730 goto x970;
2731 }
2732 /* ===============================================================
2733 Move the pivot element to the front of its row
2734 and to the top of its column.
2735 =============================================================== */
2736 LPIVR = LUSOL->locr[IBEST];
2737 LPIVR1 = LPIVR+1;
2738 LPIVR2 = LPIVR+NELIM;
2739 for(L = LPIVR; L <= LPIVR2; L++) {
2740 if(LUSOL->indr[L]==JBEST)
2741 break;
2742 }
2743
2744 LUSOL->indr[L] = LUSOL->indr[LPIVR];
2745 LUSOL->indr[LPIVR] = JBEST;
2746 LPIVC = LUSOL->locc[JBEST];
2747 LPIVC1 = LPIVC+1;
2748 LPIVC2 = LPIVC+MELIM;
2749 for(L = LPIVC; L <= LPIVC2; L++) {
2750 if(LUSOL->indc[L]==IBEST)
2751 break;
2752 }
2753 LUSOL->indc[L] = LUSOL->indc[LPIVC];
2754 LUSOL->indc[LPIVC] = IBEST;
2755 ABEST = LUSOL->a[L];
2756 LUSOL->a[L] = LUSOL->a[LPIVC];
2757 LUSOL->a[LPIVC] = ABEST;
2758 if(!KEEPLU)
2759 /* Store just the diagonal of U, in natural order.
2760 !! a[ldiagU + nrowu] = abest ! This was in pivot order. */
2761 LUSOL->diagU[JBEST] = ABEST;
2762
2763 /* ==============================================================
2764 Delete pivot col from heap.
2765 Hk tells us where it is in the heap.
2766 ============================================================== */
2767 if(TCP) {
2768 #ifdef ClassicHamaxR
2769 KBEST = HK[JBEST];
2770 HDELETE(HA,HJ,HK,&HLEN,KBEST,&H);
2771 #else
2772 KBEST = LUSOL->Hk[JBEST];
2773 HDELETE(LUSOL->Ha,LUSOL->Hj,LUSOL->Hk,&HLEN,KBEST,&H);
2774 #endif
2775 HOPS += H;
2776 }
2777 /* ===============================================================
2778 Delete the pivot row from the column file
2779 and store it as the next row of U.
2780 set indr(lu) = 0 to initialize jfill ptrs on columns of D,
2781 indc(lu) = lenj to save the original column lengths.
2782 =============================================================== */
2783 LUSOL->a[LU1] = ABEST;
2784 LUSOL->indr[LU1] = JBEST;
2785 LUSOL->indc[LU1] = NROWD;
2786 LU = LU1;
2787 DIAG = fabs(ABEST);
2788 SETMAX(*UMAX,DIAG);
2789 SETMAX(*DUMAX,DIAG);
2790 SETMIN(*DUMIN,DIAG);
2791 for(LR = LPIVR1; LR <= LPIVR2; LR++) {
2792 LU++;
2793 J = LUSOL->indr[LR];
2794 LENJ = LUSOL->lenc[J];
2795 LUSOL->lenc[J] = LENJ-1;
2796 LC1 = LUSOL->locc[J];
2797 LAST = LC1+LUSOL->lenc[J];
2798 for(L = LC1; L <= LAST; L++) {
2799 if(LUSOL->indc[L]==IBEST)
2800 break;
2801 }
2802 LUSOL->a[LU] = LUSOL->a[L];
2803 LUSOL->indr[LU] = 0;
2804 LUSOL->indc[LU] = LENJ;
2805 SETMAX(*UMAX,fabs(LUSOL->a[LU]));
2806 LUSOL->a[L] = LUSOL->a[LAST];
2807 LUSOL->indc[L] = LUSOL->indc[LAST];
2808 /* Free entry */
2809 LUSOL->indc[LAST] = 0;
2810 /* ??? if (j .eq. jlast) lcol = lcol - 1 */
2811 }
2812 /* ===============================================================
2813 Delete the pivot column from the row file
2814 and store the nonzeros of the next column of L.
2815 Set indc(ll) = 0 to initialize markl(*) markers,
2816 indr(ll) = 0 to initialize ifill(*) row fill-in cntrs,
2817 indc(ls) = leni to save the original row lengths,
2818 indr(ls) = iqloc(i) to save parts of iqloc(*),
2819 iqloc(i) = lsave - ls to point to the nonzeros of L
2820 = -1, -2, -3, ... in mark(*).
2821 =============================================================== */
2822 LUSOL->indc[LSAVE] = NCOLD;
2823 if(MELIM==0)
2824 goto x700;
2825 LL = LL1-1;
2826 LS = LSAVE;
2827 ABEST = ONE/ABEST;
2828 for(LC = LPIVC1; LC <= LPIVC2; LC++) {
2829 LL++;
2830 LS++;
2831 I = LUSOL->indc[LC];
2832 LENI = LUSOL->lenr[I];
2833 LUSOL->lenr[I] = LENI-1;
2834 LR1 = LUSOL->locr[I];
2835 LAST = LR1+LUSOL->lenr[I];
2836 for(L = LR1; L <= LAST; L++) {
2837 if(LUSOL->indr[L]==JBEST)
2838 break;
2839 }
2840 LUSOL->indr[L] = LUSOL->indr[LAST];
2841 /* Free entry */
2842 LUSOL->indr[LAST] = 0;
2843 LUSOL->a[LL] = -LUSOL->a[LC]*ABEST;
2844 LIJ = fabs(LUSOL->a[LL]);
2845 SETMAX(*LMAX,LIJ);
2846 LUSOL->indc[LL] = 0;
2847 LUSOL->indr[LL] = 0;
2848 LUSOL->indc[LS] = LENI;
2849 LUSOL->indr[LS] = LUSOL->iqloc[I];
2850 LUSOL->iqloc[I] = LSAVE-LS;
2851 }
2852 /* ===============================================================
2853 Do the Gaussian elimination.
2854 This involves adding a multiple of the pivot column
2855 to all other columns in the pivot row.
2856 Sometimes more than one call to lu1gau is needed to allow
2857 compression of the column file.
2858 lfirst says which column the elimination should start with.
2859 minfre is a bound on the storage needed for any one column.
2860 lu points to off-diagonals of u.
2861 nfill keeps track of pending fill-in in the row file.
2862 =============================================================== */
2863 if(NELIM==0)
2864 goto x700;
2865 LFIRST = LPIVR1;
2866 MINFRE = MLEFT+NSPARE;
2867 LU = 1;
2868 NFILL = 0;
2869
2870 x400:
2871 #ifdef UseTimer
2872 timer ( "start", eltime );
2873 #endif
2874 LU1GAU(LUSOL, MELIM,NSPARE,SMALL,LPIVC1,LPIVC2,&LFIRST,LPIVR2,
2875 LFREE,MINFRE,ILAST,&JLAST,&LROW,&LCOL,&LU,&NFILL,
2876 LUSOL->iqloc, LUSOL->a+LL1-LUSOL_ARRAYOFFSET,
2877 LUSOL->indc+LL1-LUSOL_ARRAYOFFSET, LUSOL->a+LU1-LUSOL_ARRAYOFFSET,
2878 LUSOL->indr+LL1-LUSOL_ARRAYOFFSET, LUSOL->indr+LU1-LUSOL_ARRAYOFFSET);
2879 #ifdef UseTimer
2880 timer ( "finish", eltime );
2881 #endif
2882 if(LFIRST>0) {
2883 /* The elimination was interrupted.
2884 Compress the column file and try again.
2885 lfirst, lu and nfill have appropriate new values. */
2886 LU1REC(LUSOL, LUSOL->n,TRUE,&LCOL,
2887 LUSOL->indc,LUSOL->lenc,LUSOL->locc);
2888 LFILE = LCOL;
2889 JLAST = LUSOL->indc[LCOL+1];
2890 LPIVC = LUSOL->locc[JBEST];
2891 LPIVC1 = LPIVC+1;
2892 LPIVC2 = LPIVC+MELIM;
2893 NFREE = LFREE-LCOL;
2894 if(NFREE<MINFRE) {
2895 goto x970;
2896 }
2897 goto x400;
2898 }
2899 /* ===============================================================
2900 The column file has been fully updated.
2901 Deal with any pending fill-in in the row file.
2902 =============================================================== */
2903 if(NFILL>0) {
2904 /* Compress the row file if necessary.
2905 lu1gau has set nfill to be the number of pending fill-ins
2906 plus the current length of any rows that need to be moved. */
2907 MINFRE = NFILL;
2908 NFREE = LFREE-LROW;
2909 if(NFREE<MINFRE) {
2910 LU1REC(LUSOL, LUSOL->m,FALSE,&LROW,
2911 LUSOL->indr,LUSOL->lenr,LUSOL->locr);
2912 LFILE = LROW;
2913 ILAST = LUSOL->indr[LROW+1];
2914 LPIVR = LUSOL->locr[IBEST];
2915 LPIVR1 = LPIVR+1;
2916 LPIVR2 = LPIVR+NELIM;
2917 NFREE = LFREE-LROW;
2918 if(NFREE<MINFRE) {
2919 goto x970;
2920 }
2921 }
2922 /* Move rows that have pending fill-in to end of the row file.
2923 Then insert the fill-in. */
2924 LU1PEN(LUSOL, NSPARE,&ILAST,
2925 LPIVC1,LPIVC2,LPIVR1,LPIVR2,
2926 &LROW,LUSOL->indr+LL1-LUSOL_ARRAYOFFSET,LUSOL->indr+LU1-LUSOL_ARRAYOFFSET);
2927 }
2928 /* ===============================================================
2929 Restore the saved values of iqloc.
2930 Insert the correct indices for the col of L and the row of U.
2931 =============================================================== */
2932 x700:
2933 LUSOL->lenr[IBEST] = 0;
2934 LUSOL->lenc[JBEST] = 0;
2935 LL = LL1-1;
2936 LS = LSAVE;
2937 for(LC = LPIVC1; LC <= LPIVC2; LC++) {
2938 LL++;
2939 LS++;
2940 I = LUSOL->indc[LC];
2941 LUSOL->iqloc[I] = LUSOL->indr[LS];
2942 LUSOL->indc[LL] = I;
2943 LUSOL->indr[LL] = IBEST;
2944 }
2945 LU = LU1-1;
2946 for(LR = LPIVR; LR <= LPIVR2; LR++) {
2947 LU++;
2948 LUSOL->indr[LU] = LUSOL->indr[LR];
2949 }
2950 /* ===============================================================
2951 Free the space occupied by the pivot row
2952 and update the column permutation.
2953 Then free the space occupied by the pivot column
2954 and update the row permutation.
2955 nzchng is found in both calls to lu1pq2, but we use it only
2956 after the second.
2957 =============================================================== */
2958 LU1PQ2(LUSOL, NCOLD, &NZCHNG,
2959 LUSOL->indr+LPIVR-LUSOL_ARRAYOFFSET,
2960 LUSOL->indc+LU1-LUSOL_ARRAYOFFSET, LUSOL->lenc,
2961 LUSOL->iqloc, LUSOL->iq, LUSOL->iqinv);
2962 LU1PQ2(LUSOL, NROWD, &NZCHNG,
2963 LUSOL->indc+LPIVC-LUSOL_ARRAYOFFSET,
2964 LUSOL->indc+LSAVE-LUSOL_ARRAYOFFSET, LUSOL->lenr,
2965 LUSOL->iploc, LUSOL->ip, LUSOL->ipinv);
2966 NZLEFT += NZCHNG;
2967
2968 /* ===============================================================
2969 lu1mxr resets Amaxr(i) in each modified row i.
2970 lu1mxc moves the largest aij to the top of each modified col j.
2971 28 Jun 2002: Note that cols of L have an implicit diag of 1.0,
2972 so lu1mxr is called with ll1, not ll1+1, whereas
2973 lu1mxc is called with lu1+1.
2974 =============================================================== */
2975 if(UTRI && TPP) {
2976 /* Relax -- we're not keeping big elements at the top yet. */
2977 }
2978 else {
2979 if(TRP && MELIM>0)
2980 #ifdef ClassicHamaxR
2981 LU1MXR(LUSOL, LL1,LL,LUSOL->indc,AMAXR);
2982 #else
2983 LU1MXR(LUSOL, LL1,LL,LUSOL->indc,LUSOL->amaxr);
2984 #endif
2985
2986 if(NELIM>0) {
2987 LU1MXC(LUSOL, LU1+1,LU,LUSOL->indr);
2988 /* Update modified columns in heap */
2989 if(TCP) {
2990 for(KK = LU1+1; KK <= LU; KK++) {
2991 J = LUSOL->indr[KK];
2992 #ifdef ClassicHamaxR
2993 K = HK[J];
2994 #else
2995 K = LUSOL->Hk[J];
2996 #endif
2997 /* Biggest aij in column j */
2998 V = fabs(LUSOL->a[LUSOL->locc[J]]);
2999 #ifdef ClassicHamaxR
3000 HCHANGE(HA,HJ,HK,HLEN,K,V,J,&H);
3001 #else
3002 HCHANGE(LUSOL->Ha,LUSOL->Hj,LUSOL->Hk,HLEN,K,V,J,&H);
3003 #endif
3004 HOPS += H;
3005 }
3006 }
3007 }
3008 }
3009 /* ===============================================================
3010 Negate lengths of pivot row and column so they will be
3011 eliminated during compressions.
3012 =============================================================== */
3013 LUSOL->lenr[IBEST] = -NCOLD;
3014 LUSOL->lenc[JBEST] = -NROWD;
3015
3016 /* Test for fatal bug: row or column lists overwriting L and U. */
3017 if(LROW>LSAVE || LCOL>LSAVE)
3018 goto x980;
3019
3020 /* Reset the file lengths if pivot row or col was at the end. */
3021 if(IBEST==ILAST)
3022 LROW = LUSOL->locr[IBEST];
3023
3024 if(JBEST==JLAST)
3025 LCOL = LUSOL->locc[JBEST];
3026
3027 }
3028 /* ------------------------------------------------------------------
3029 End of main loop.
3030 ------------------------------------------------------------------
3031 ------------------------------------------------------------------
3032 Normal exit.
3033 Move empty rows and cols to the end of ip, iq.
3034 Then finish with a dense LU if necessary.
3035 ------------------------------------------------------------------ */
3036 x900:
3037 *INFORM = LUSOL_INFORM_LUSUCCESS;
3038 LU1PQ3(LUSOL, LUSOL->m,LUSOL->lenr,LUSOL->ip,LUSOL->ipinv,&MRANK);
3039 LU1PQ3(LUSOL, LUSOL->n,LUSOL->lenc,LUSOL->iq,LUSOL->iqinv,NRANK);
3040 SETMIN(*NRANK, MRANK);
3041 if(DENSLU) {
3042 #ifdef UseTimer
3043 timer ( "start", 17 );
3044 #endif
3045 LU1FUL(LUSOL, LEND,LU1,TPP,MLEFT,NLEFT,*NRANK,NROWU,LENL,LENU,
3046 &NSING,KEEPLU,SMALL,LUSOL->a+LD-LUSOL_ARRAYOFFSET,LUSOL->locr);
3047 /* *** 21 Dec 1994: Bug in next line.
3048 *** nrank = nrank - nsing */
3049 *NRANK = MINMN-NSING;
3050 #ifdef UseTimer
3051 timer ( "finish", 17 );
3052 #endif
3053 }
3054 *MINLEN = (*LENL)+(*LENU)+2*(LUSOL->m+LUSOL->n);
3055 goto x990;
3056 /* Not enough space free after a compress.
3057 Set minlen to an estimate of the necessary value of lena. */
3058 x970:
3059 *INFORM = LUSOL_INFORM_ANEEDMEM;
3060 *MINLEN = LENA2+LFILE+2*(LUSOL->m+LUSOL->n);
3061 goto x990;
3062 /* Fatal error. This will never happen!
3063 (Famous last words.) */
3064 x980:
3065 *INFORM = LUSOL_INFORM_FATALERR;
3066 goto x990;
3067 /* Fatal error with TSP. Diagonal pivot not found. */
3068 x985:
3069 *INFORM = LUSOL_INFORM_NOPIVOT;
3070 /* Exit. */
3071 x990:
3072 #ifdef UseTimer
3073 timer ( "finish", 3 );
3074 #endif
3075 ;
3076 }
3077
3078
3079 /* ==================================================================
3080 lu1fac computes a factorization A = L*U, where A is a sparse
3081 matrix with m rows and n columns, P*L*P' is lower triangular
3082 and P*U*Q is upper triangular for certain permutations P, Q
3083 (which are returned in the arrays ip, iq).
3084 Stability is ensured by limiting the size of the elements of L.
3085 The nonzeros of A are input via the parallel arrays a, indc, indr,
3086 which should contain nelem entries of the form aij, i, j
3087 in any order. There should be no duplicate pairs i, j.
3088
3089 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3090 + Beware !!! The row indices i must be in indc, +
3091 + and the column indices j must be in indr. +
3092 + (Not the other way round!) +
3093 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3094
3095 It does not matter if some of the entries in a(*) are zero.
3096 Entries satisfying abs( a(i) ) .le. parmlu(3) are ignored.
3097 Other parameters in luparm and parmlu are described below.
3098 The matrix A may be singular. On exit, nsing = luparm(11) gives
3099 the number of apparent singularities. This is the number of
3100 "small" diagonals of the permuted factor U, as judged by
3101 the input tolerances Utol1 = parmlu(4) and Utol2 = parmlu(5).
3102 The diagonal element diagj associated with column j of A is
3103 "small" if
3104 abs( diagj ) .le. Utol1
3105 or
3106 abs( diagj ) .le. Utol2 * max( uj ),
3107 where max( uj ) is the maximum element in the j-th column of U.
3108 The position of such elements is returned in w(*). In general,
3109 w(j) = + max( uj ), but if column j is a singularity,
3110 w(j) = - max( uj ). Thus, w(j) .le. 0 if column j appears to be
3111 dependent on the other columns of A.
3112 NOTE: lu1fac (like certain other sparse LU packages) does not
3113 treat dense columns efficiently. This means it will be slow
3114 on "arrow matrices" of the form
3115 A = (x a)
3116 ( x b)
3117 ( x c)
3118 ( x d)
3119 (x x x x e)
3120 if the numerical values in the dense column allow it to be
3121 chosen LATE in the pivot order.
3122 With TPP (Threshold Partial Pivoting), the dense column is
3123 likely to be chosen late.
3124 With TCP (Threshold Complete Pivoting), if any of a,b,c,d
3125 is significantly larger than other elements of A, it will
3126 be chosen as the first pivot and the dense column will be
3127 eliminated, giving reasonably sparse factors.
3128 However, if element e is so big that TCP chooses it, the factors
3129 will become dense. (It's hard to win on these examples!)
3130 ------------------------------------------------------------------
3131
3132 Notes on the array names
3133 ------------------------
3134 During the LU factorization, the sparsity pattern of the matrix
3135 being factored is stored twice: in a column list and a row list.
3136 The column list is ( a, indc, locc, lenc )
3137 where
3138 a(*) holds the nonzeros,
3139 indc(*) holds the indices for the column list,
3140 locc(j) points to the start of column j in a(*) and indc(*),
3141 lenc(j) is the number of nonzeros in column j.
3142 The row list is ( indr, locr, lenr )
3143 where
3144 indr(*) holds the indices for the row list,
3145 locr(i) points to the start of row i in indr(*),
3146 lenr(i) is the number of nonzeros in row i.
3147 At all stages of the LU factorization, ip contains a complete
3148 row permutation. At the start of stage k, ip(1), ..., ip(k-1)
3149 are the first k-1 rows of the final row permutation P.
3150 The remaining rows are stored in an ordered list
3151 ( ip, iploc, ipinv )
3152 where
3153 iploc(nz) points to the start in ip(*) of the set of rows
3154 that currently contain nz nonzeros,
3155 ipinv(i) points to the position of row i in ip(*).
3156 For example,
3157 iploc(1) = k (and this is where rows of length 1 {),
3158 iploc(2) = k+p if there are p rows of length 1
3159 (and this is where rows of length 2 {).
3160 Similarly for iq, iqloc, iqinv.
3161 ---------------------------------------------------------------------
3162 INPUT PARAMETERS
3163 m (not altered) is the number of rows in A.
3164 n (not altered) is the number of columns in A.
3165 nelem (not altered) is the number of matrix entries given in
3166 the arrays a, indc, indr.
3167 lena (not altered) is the dimension of a, indc, indr.
3168 This should be significantly larger than nelem.
3169 Typically one should have
3170 lena > max( 2*nelem, 10*m, 10*n, 10000 )
3171 but some applications may need more.
3172 On machines with virtual memory it is safe to have
3173 lena "far bigger than necessary", since not all of the
3174 arrays will be used.
3175 a (overwritten) contains entries Aij in a(1:nelem).
3176 indc (overwritten) contains the indices i in indc(1:nelem).
3177 indr (overwritten) contains the indices j in indr(1:nelem).
3178 luparm input parameters: Typical value
3179 luparm( 1) = nout File number for printed messages. 6
3180 luparm( 2) = lprint Print level. 0
3181 < 0 suppresses output.
3182 = 0 gives error messages.
3183 >= 10 gives statistics about the LU factors.
3184 >= 50 gives debug output from lu1fac
3185 (the pivot row and column and the
3186 no. of rows and columns involved at
3187 each elimination step).
3188 luparm( 3) = maxcol lu1fac: maximum number of columns 5
3189 searched allowed in a Markowitz-type
3190 search for the next pivot element.
3191 For some of the factorization, the
3192 number of rows searched is
3193 maxrow = maxcol - 1.
3194 luparm( 6) = 0 => TPP: Threshold Partial Pivoting. 0
3195 = 1 => TRP: Threshold Rook Pivoting.
3196 = 2 => TCP: Threshold Complete Pivoting.
3197 = 3 => TSP: Threshold Symmetric Pivoting.
3198 = 4 => TDP: Threshold Diagonal Pivoting.
3199 (TDP not yet implemented).
3200 TRP and TCP are more expensive than TPP but
3201 more stable and better at revealing rank.
3202 Take care with setting parmlu(1), especially
3203 with TCP.
3204 NOTE: TSP and TDP are for symmetric matrices
3205 that are either definite or quasi-definite.
3206 TSP is effectively TRP for symmetric matrices.
3207 TDP is effectively TCP for symmetric matrices.
3208 luparm( 8) = keepLU lu1fac: keepLU = 1 means the numerical 1
3209 factors will be computed if possible.
3210 keepLU = 0 means L and U will be discarded
3211 but other information such as the row and
3212 column permutations will be returned.
3213 The latter option requires less storage.
3214 parmlu input parameters: Typical value
3215 parmlu( 1) = Ltol1 Max Lij allowed during Factor.
3216 TPP 10.0 or 100.0
3217 TRP 4.0 or 10.0
3218 TCP 5.0 or 10.0
3219 TSP 4.0 or 10.0
3220 With TRP and TCP (Rook and Complete Pivoting),
3221 values less than 25.0 may be expensive
3222 on badly scaled data. However,
3223 values less than 10.0 may be needed
3224 to obtain a reliable rank-revealing
3225 factorization.
3226 parmlu( 2) = Ltol2 Max Lij allowed during Updates. 10.0
3227 during updates.
3228 parmlu( 3) = small Absolute tolerance for eps**0.8 = 3.0d-13
3229 treating reals as zero.
3230 parmlu( 4) = Utol1 Absolute tol for flagging eps**0.67= 3.7d-11
3231 small diagonals of U.
3232 parmlu( 5) = Utol2 Relative tol for flagging eps**0.67= 3.7d-11
3233 small diagonals of U.
3234 (eps = machine precision)
3235 parmlu( 6) = Uspace Factor limiting waste space in U. 3.0
3236 In lu1fac, the row or column lists
3237 are compressed if their length
3238 exceeds Uspace times the length of
3239 either file after the last compression.
3240 parmlu( 7) = dens1 The density at which the Markowitz 0.3
3241 pivot strategy should search maxcol
3242 columns and no rows.
3243 (Use 0.3 unless you are experimenting
3244 with the pivot strategy.)
3245 parmlu( 8) = dens2 the density at which the Markowitz 0.5
3246 strategy should search only 1 column,
3247 or (if storage is available)
3248 the density at which all remaining
3249 rows and columns will be processed
3250 by a dense LU code.
3251 For example, if dens2 = 0.1 and lena is
3252 large enough, a dense LU will be used
3253 once more than 10 per cent of the
3254 remaining matrix is nonzero.
3255
3256 OUTPUT PARAMETERS
3257 a, indc, indr contain the nonzero entries in the LU factors of A.
3258 If keepLU = 1, they are in a form suitable for use
3259 by other parts of the LUSOL package, such as lu6sol.
3260 U is stored by rows at the start of a, indr.
3261 L is stored by cols at the end of a, indc.
3262 If keepLU = 0, only the diagonals of U are stored, at the
3263 end of a.
3264 ip, iq are the row and column permutations defining the
3265 pivot order. For example, row ip(1) and column iq(1)
3266 defines the first diagonal of U.
3267 lenc(1:numl0) contains the number of entries in nontrivial
3268 columns of L (in pivot order).
3269 lenr(1:m) contains the number of entries in each row of U
3270 (in original order).
3271 locc(1:n) = 0 (ready for the LU update routines).
3272 locr(1:m) points to the beginning of the rows of U in a, indr.
3273 iploc, iqloc, ipinv, iqinv are undefined.
3274 w indicates singularity as described above.
3275 inform = 0 if the LU factors were obtained successfully.
3276 = 1 if U appears to be singular, as judged by lu6chk.
3277 = 3 if some index pair indc(l), indr(l) lies outside
3278 the matrix dimensions 1:m , 1:n.
3279 = 4 if some index pair indc(l), indr(l) duplicates
3280 another such pair.
3281 = 7 if the arrays a, indc, indr were not large enough.
3282 Their length "lena" should be increase to at least
3283 the value "minlen" given in luparm(13).
3284 = 8 if there was some other fatal error. (Shouldn't happen!)
3285 = 9 if no diagonal pivot could be found with TSP or TDP.
3286 The matrix must not be sufficiently definite
3287 or quasi-definite.
3288 luparm output parameters:
3289 luparm(10) = inform Return code from last call to any LU routine.
3290 luparm(11) = nsing No. of singularities marked in the
3291 output array w(*).
3292 luparm(12) = jsing Column index of last singularity.
3293 luparm(13) = minlen Minimum recommended value for lena.
3294 luparm(14) = maxlen ?
3295 luparm(15) = nupdat No. of updates performed by the lu8 routines.
3296 luparm(16) = nrank No. of nonempty rows of U.
3297 luparm(17) = ndens1 No. of columns remaining when the density of
3298 the matrix being factorized reached dens1.
3299 luparm(18) = ndens2 No. of columns remaining when the density of
3300 the matrix being factorized reached dens2.
3301 luparm(19) = jumin The column index associated with DUmin.
3302 luparm(20) = numL0 No. of columns in initial L.
3303 luparm(21) = lenL0 Size of initial L (no. of nonzeros).
3304 luparm(22) = lenU0 Size of initial U.
3305 luparm(23) = lenL Size of current L.
3306 luparm(24) = lenU Size of current U.
3307 luparm(25) = lrow Length of row file.
3308 luparm(26) = ncp No. of compressions of LU data structures.
3309 luparm(27) = mersum lu1fac: sum of Markowitz merit counts.
3310 luparm(28) = nUtri lu1fac: triangular rows in U.
3311 luparm(29) = nLtri lu1fac: triangular rows in L.
3312 luparm(30) =
3313 parmlu output parameters:
3314 parmlu(10) = Amax Maximum element in A.
3315 parmlu(11) = Lmax Maximum multiplier in current L.
3316 parmlu(12) = Umax Maximum element in current U.
3317 parmlu(13) = DUmax Maximum diagonal in U.
3318 parmlu(14) = DUmin Minimum diagonal in U.
3319 parmlu(15) = Akmax Maximum element generated at any stage
3320 during TCP factorization.
3321 parmlu(16) = growth TPP: Umax/Amax TRP, TCP, TSP: Akmax/Amax
3322 parmlu(17) =
3323 parmlu(18) =
3324 parmlu(19) =
3325 parmlu(20) = resid lu6sol: residual after solve with U or U'.
3326 ...
3327 parmlu(30) =
3328 ------------------------------------------------------------------
3329 00 Jun 1983 Original version.
3330 00 Jul 1987 nrank saved in luparm(16).
3331 12 Apr 1989 ipinv, iqinv added as workspace.
3332 26 Apr 1989 maxtie replaced by maxcol in Markowitz search.
3333 16 Mar 1992 jumin saved in luparm(19).
3334 10 Jun 1992 lu1fad has to move empty rows and cols to the bottom
3335 (via lu1pq3) before doing the dense LU.
3336 12 Jun 1992 Deleted dense LU (lu1ful, lu1vlu).
3337 25 Oct 1993 keepLU implemented.
3338 07 Feb 1994 Added new dense LU (lu1ful, lu1den).
3339 21 Dec 1994 Bugs fixed in lu1fad (nrank) and lu1ful (ipvt).
3340 08 Aug 1995 Use ip instead of w as parameter to lu1or3 (for F90).
3341 13 Sep 2000 TPP and TCP options implemented.
3342 17 Oct 2000 Fixed troubles due to A = empty matrix (Todd Munson).
3343 01 Dec 2000 Save Lmax, Umax, etc. after both lu1fad and lu6chk.
3344 lu1fad sets them when keepLU = false.
3345 lu6chk sets them otherwise, and includes items
3346 from the dense LU.
3347 11 Mar 2001 lu6chk now looks at diag(U) when keepLU = false.
3348 26 Apr 2002 New TCP implementation using heap routines to
3349 store largest element in each column.
3350 New workspace arrays Ha, Hj, Hk required.
3351 For compatibility, borrow space from a, indc, indr
3352 rather than adding new input parameters.
3353 01 May 2002 lu1den changed to lu1DPP (dense partial pivoting).
3354 lu1DCP implemented (dense complete pivoting).
3355 Both TPP and TCP now switch to dense mode and end.
3356 ================================================================== */
3357 void LU1FAC(LUSOLrec *LUSOL, int *INFORM)
3358 {
3359 MYBOOL KEEPLU, TCP, TPP, TRP, TSP;
3360 int LPIV, NELEM0, LPRINT, MINLEN, NUML0, LENL, LENU, LROW, MERSUM,
3361 NUTRI, NLTRI, NDENS1, NDENS2, NRANK, NSING, JSING, JUMIN, NUMNZ, LERR,
3362 LU, LL, LM, LTOPL, K, I, LENUK, J, LENLK, IDUMMY, LLSAVE, NMOVE, L2, L, NCP, NBUMP;
3363 #ifdef ClassicHamaxR
3364 int LENH, LENA2, LOCH, LMAXR;
3365 #endif
3366
3367 REAL LMAX, LTOL, SMALL, AMAX, UMAX, DUMAX, DUMIN, AKMAX, DM, DN, DELEM, DENSTY,
3368 AGRWTH, UGRWTH, GROWTH, CONDU, DINCR, AVGMER;
3369
3370 /* Free row-based version of L0 (regenerated by LUSOL_btran). */
3371 if(LUSOL->L0 != NULL)
3372 LUSOL_matfree(&(LUSOL->L0));
3373
3374 /* Grab relevant input parameters. */
3375 NELEM0 = LUSOL->nelem;
3376 LPRINT = LUSOL->luparm[LUSOL_IP_PRINTLEVEL];
3377 LPIV = LUSOL->luparm[LUSOL_IP_PIVOTTYPE];
3378 KEEPLU = (MYBOOL) (LUSOL->luparm[LUSOL_IP_KEEPLU]!=FALSE);
3379 /* Limit on size of Lij */
3380 LTOL = LUSOL->parmlu[LUSOL_RP_FACTORMAX_Lij];
3381 /* Drop tolerance */
3382 SMALL = LUSOL->parmlu[LUSOL_RP_ZEROTOLERANCE];
3383 TPP = (MYBOOL) (LPIV==LUSOL_PIVMOD_TPP);
3384 TRP = (MYBOOL) (LPIV==LUSOL_PIVMOD_TRP);
3385 TCP = (MYBOOL) (LPIV==LUSOL_PIVMOD_TCP);
3386 TSP = (MYBOOL) (LPIV==LUSOL_PIVMOD_TSP);
3387 /* Initialize output parameters. */
3388 *INFORM = LUSOL_INFORM_LUSUCCESS;
3389 LERR = 0;
3390 MINLEN = LUSOL->nelem + 2*(LUSOL->m+LUSOL->n);
3391 NUML0 = 0;
3392 LENL = 0;
3393 LENU = 0;
3394 LROW = 0;
3395 MERSUM = 0;
3396 NUTRI = LUSOL->m;
3397 NLTRI = 0;
3398 NDENS1 = 0;
3399 NDENS2 = 0;
3400 NRANK = 0;
3401 NSING = 0;
3402 JSING = 0;
3403 JUMIN = 0;
3404 AMAX = ZERO;
3405 LMAX = ZERO;
3406 UMAX = ZERO;
3407 DUMAX = ZERO;
3408 DUMIN = ZERO;
3409 AKMAX = ZERO;
3410
3411 /* Float version of dimensions. */
3412 DM = LUSOL->m;
3413 DN = LUSOL->n;
3414 DELEM = LUSOL->nelem;
3415
3416 /* Initialize workspace parameters. */
3417 LUSOL->luparm[LUSOL_IP_COMPRESSIONS_LU] = 0;
3418 if(LUSOL->lena < MINLEN) {
3419 if(!LUSOL_realloc_a(LUSOL, MINLEN))
3420 goto x970;
3421 }
3422
3423 /* ------------------------------------------------------------------
3424 Organize the aij's in a, indc, indr.
3425 lu1or1 deletes small entries, tests for illegal i,j's,
3426 and counts the nonzeros in each row and column.
3427 lu1or2 reorders the elements of A by columns.
3428 lu1or3 uses the column list to test for duplicate entries
3429 (same indices i,j).
3430 lu1or4 constructs a row list from the column list.
3431 ------------------------------------------------------------------ */
3432 LU1OR1(LUSOL, SMALL,&AMAX,&NUMNZ,&LERR,INFORM);
3433 if(LPRINT>=LUSOL_MSG_STATISTICS) {
3434 DENSTY = (100*DELEM)/(DM*DN);
3435 LUSOL_report(LUSOL, 0, "m:%6d %c n:%6d nzcount:%9d Amax:%g Density:%g\n",
3436 LUSOL->m, relationChar(LUSOL->m, LUSOL->n), LUSOL->n,
3437 LUSOL->nelem, AMAX, DENSTY);
3438 }
3439 if(*INFORM!=LUSOL_INFORM_LUSUCCESS)
3440 goto x930;
3441 LUSOL->nelem = NUMNZ;
3442 LU1OR2(LUSOL);
3443 LU1OR3(LUSOL, &LERR,INFORM);
3444 if(*INFORM!=LUSOL_INFORM_LUSUCCESS)
3445 goto x940;
3446 LU1OR4(LUSOL);
3447 /* ------------------------------------------------------------------
3448 Set up lists of rows and columns with equal numbers of nonzeros,
3449 using indc(*) as workspace.
3450 ------------------------------------------------------------------ */
3451 LU1PQ1(LUSOL, LUSOL->m,LUSOL->n,LUSOL->lenr,
3452 LUSOL->ip,LUSOL->iploc,LUSOL->ipinv,
3453 LUSOL->indc+LUSOL->nelem); /* LUSOL_ARRAYOFFSET implied */
3454 LU1PQ1(LUSOL, LUSOL->n,LUSOL->m,LUSOL->lenc,
3455 LUSOL->iq,LUSOL->iqloc,LUSOL->iqinv,
3456 LUSOL->indc+LUSOL->nelem); /* LUSOL_ARRAYOFFSET implied */
3457 /* ------------------------------------------------------------------
3458 For TCP, Ha, Hj, Hk are allocated separately, similarly amaxr
3459 for TRP. Then compute the factorization A = L*U.
3460 ------------------------------------------------------------------ */
3461 #ifdef ClassicHamaxR
3462 if(TPP || TSP) {
3463 LENH = 1;
3464 LENA2 = LUSOL->lena;
3465 LOCH = LUSOL->lena;
3466 LMAXR = 1;
3467 }
3468 else if(TRP) {
3469 LENH = 1; /* Dummy */
3470 LENA2 = LUSOL->lena-LUSOL->m; /* Reduced length of a */
3471 LOCH = LUSOL->lena; /* Dummy */
3472 LMAXR = LENA2+1; /* Start of Amaxr in a */
3473 }
3474 else if(TCP) {
3475 LENH = LUSOL->n; /* Length of heap */
3476 LENA2 = LUSOL->lena-LENH; /* Reduced length of a, indc, indr */
3477 LOCH = LENA2+1; /* Start of Ha, Hj, Hk in a, indc, indr */
3478 LMAXR = 1; /* Dummy */
3479 }
3480 LU1FAD(LUSOL,
3481 LENA2,LENH,
3482 LUSOL->a+LOCH-LUSOL_ARRAYOFFSET,
3483 LUSOL->indc+LOCH-LUSOL_ARRAYOFFSET,
3484 LUSOL->indr+LOCH-LUSOL_ARRAYOFFSET,
3485 LUSOL->a+LMAXR-LUSOL_ARRAYOFFSET,
3486 INFORM,&LENL,&LENU,
3487 &MINLEN,&MERSUM,&NUTRI,&NLTRI,&NDENS1,&NDENS2,
3488 &NRANK,&LMAX,&UMAX,&DUMAX,&DUMIN,&AKMAX);
3489 #else
3490 LU1FAD(LUSOL,
3491 INFORM,&LENL,&LENU,
3492 &MINLEN,&MERSUM,&NUTRI,&NLTRI,&NDENS1,&NDENS2,
3493 &NRANK,&LMAX,&UMAX,&DUMAX,&DUMIN,&AKMAX);
3494 #endif
3495 LUSOL->luparm[LUSOL_IP_RANK_U] = NRANK;
3496 LUSOL->luparm[LUSOL_IP_NONZEROS_L] = LENL;
3497 if(*INFORM==LUSOL_INFORM_ANEEDMEM)
3498 goto x970;
3499 if(*INFORM==LUSOL_INFORM_NOPIVOT)
3500 goto x985;
3501 if(*INFORM>LUSOL_INFORM_LUSUCCESS)
3502 goto x980;
3503 if(KEEPLU) {
3504 /* ---------------------------------------------------------------
3505 The LU factors are at the top of a, indc, indr,
3506 with the columns of L and the rows of U in the order
3507 ( free ) ... ( u3 ) ( l3 ) ( u2 ) ( l2 ) ( u1 ) ( l1 ).
3508 Starting with ( l1 ) and ( u1 ), move the rows of U to the
3509 left and the columns of L to the right, giving
3510 ( u1 ) ( u2 ) ( u3 ) ... ( free ) ... ( l3 ) ( l2 ) ( l1 ).
3511 Also, set numl0 = the number of nonempty columns of U.
3512 --------------------------------------------------------------- */
3513 LU = 0;
3514 LL = LUSOL->lena+1;
3515 #ifdef ClassicHamaxR
3516 LM = LENA2+1;
3517 #else
3518 LM = LL;
3519 #endif
3520 LTOPL = LL-LENL-LENU;
3521 LROW = LENU;
3522 for(K = 1; K <= NRANK; K++) {
3523 I = LUSOL->ip[K];
3524 LENUK = -LUSOL->lenr[I];
3525 LUSOL->lenr[I] = LENUK;
3526 J = LUSOL->iq[K];
3527 LENLK = -LUSOL->lenc[J]-1;
3528 if(LENLK>0) {
3529 NUML0++;
3530 LUSOL->iqloc[NUML0] = LENLK;
3531 }
3532 if(LU+LENUK<LTOPL) {
3533 /* =========================================================
3534 There is room to move ( uk ). Just right-shift ( lk ).
3535 ========================================================= */
3536 for(IDUMMY = 1; IDUMMY <= LENLK; IDUMMY++) {
3537 LL--;
3538 LM--;
3539 LUSOL->a[LL] = LUSOL->a[LM];
3540 LUSOL->indc[LL] = LUSOL->indc[LM];
3541 LUSOL->indr[LL] = LUSOL->indr[LM];
3542 }
3543 }
3544 else {
3545 /* =========================================================
3546 There is no room for ( uk ) yet. We have to
3547 right-shift the whole of the remaining LU file.
3548 Note that ( lk ) ends up in the correct place.
3549 ========================================================= */
3550 LLSAVE = LL-LENLK;
3551 NMOVE = LM-LTOPL;
3552 for(IDUMMY = 1; IDUMMY <= NMOVE; IDUMMY++) {
3553 LL--;
3554 LM--;
3555 LUSOL->a[LL] = LUSOL->a[LM];
3556 LUSOL->indc[LL] = LUSOL->indc[LM];
3557 LUSOL->indr[LL] = LUSOL->indr[LM];
3558 }
3559 LTOPL = LL;
3560 LL = LLSAVE;
3561 LM = LL;
3562 }
3563 /* ======================================================
3564 Left-shift ( uk ).
3565 ====================================================== */
3566 LUSOL->locr[I] = LU+1;
3567 L2 = LM-1;
3568 LM = LM-LENUK;
3569 for(L = LM; L <= L2; L++) {
3570 LU = LU+1;
3571 LUSOL->a[LU] = LUSOL->a[L];
3572 LUSOL->indr[LU] = LUSOL->indr[L];
3573 }
3574 }
3575 /* ---------------------------------------------------------------
3576 Save the lengths of the nonempty columns of L,
3577 and initialize locc(j) for the LU update routines.
3578 --------------------------------------------------------------- */
3579 for(K = 1; K <= NUML0; K++) {
3580 LUSOL->lenc[K] = LUSOL->iqloc[K];
3581 }
3582 for(J = 1; J <= LUSOL->n; J++) {
3583 LUSOL->locc[J] = 0;
3584 }
3585 /* ---------------------------------------------------------------
3586 Test for singularity.
3587 lu6chk sets nsing, jsing, jumin, Lmax, Umax, DUmax, DUmin
3588 (including entries from the dense LU).
3589 inform = 1 if there are singularities (nsing gt 0).
3590 --------------------------------------------------------------- */
3591 LU6CHK(LUSOL, 1,LUSOL->lena,INFORM);
3592 NSING = LUSOL->luparm[LUSOL_IP_SINGULARITIES];
3593 JSING = LUSOL->luparm[LUSOL_IP_SINGULARINDEX];
3594 JUMIN = LUSOL->luparm[LUSOL_IP_COLINDEX_DUMIN];
3595 LMAX = LUSOL->parmlu[LUSOL_RP_MAXMULT_L];
3596 UMAX = LUSOL->parmlu[LUSOL_RP_MAXELEM_U];
3597 DUMAX = LUSOL->parmlu[LUSOL_RP_MAXELEM_DIAGU];
3598 DUMIN = LUSOL->parmlu[LUSOL_RP_MINELEM_DIAGU];
3599 }
3600 else {
3601 /* ---------------------------------------------------------------
3602 keepLU = 0. L and U were not kept, just the diagonals of U.
3603 lu1fac will probably be called again soon with keepLU = .true.
3604 11 Mar 2001: lu6chk revised. We can call it with keepLU = 0,
3605 but we want to keep Lmax, Umax from lu1fad.
3606 05 May 2002: Allow for TCP with new lu1DCP. Diag(U) starts
3607 below lena2, not lena. Need lena2 in next line.
3608 --------------------------------------------------------------- */
3609 #ifdef ClassicHamaxR
3610 LU6CHK(LUSOL, 1,LENA2,INFORM);
3611 #else
3612 LU6CHK(LUSOL, 1,LUSOL->lena,INFORM);
3613 #endif
3614 NSING = LUSOL->luparm[LUSOL_IP_SINGULARITIES];
3615 JSING = LUSOL->luparm[LUSOL_IP_SINGULARINDEX];
3616 JUMIN = LUSOL->luparm[LUSOL_IP_COLINDEX_DUMIN];
3617 DUMAX = LUSOL->parmlu[LUSOL_RP_MAXELEM_DIAGU];
3618 DUMIN = LUSOL->parmlu[LUSOL_RP_MINELEM_DIAGU];
3619 }
3620 goto x990;
3621 /* ------------
3622 Error exits.
3623 ------------ */
3624 x930:
3625 *INFORM = LUSOL_INFORM_ADIMERR;
3626 if(LPRINT>=LUSOL_MSG_SINGULARITY)
3627 LUSOL_report(LUSOL, 0, "lu1fac error...\nentry a[%d] has an illegal row (%d) or column (%d) index\n",
3628 LERR,LUSOL->indc[LERR],LUSOL->indr[LERR]);
3629 goto x990;
3630 x940:
3631 *INFORM = LUSOL_INFORM_ADUPLICATE;
3632 if(LPRINT>=LUSOL_MSG_SINGULARITY)
3633 LUSOL_report(LUSOL, 0, "lu1fac error...\nentry a[%d] is a duplicate with indeces indc=%d, indr=%d\n",
3634 LERR,LUSOL->indc[LERR],LUSOL->indr[LERR]);
3635 goto x990;
3636 x970:
3637 *INFORM = LUSOL_INFORM_ANEEDMEM;
3638 if(LPRINT>=LUSOL_MSG_SINGULARITY)
3639 LUSOL_report(LUSOL, 0, "lu1fac error...\ninsufficient storage; increase lena from %d to at least %d\n",
3640 LUSOL->lena, MINLEN);
3641 goto x990;
3642 x980:
3643 *INFORM = LUSOL_INFORM_FATALERR;
3644 if(LPRINT>=LUSOL_MSG_SINGULARITY)
3645 LUSOL_report(LUSOL, 0, "lu1fac error...\nfatal bug (sorry --- this should never happen)\n");
3646 goto x990;
3647 x985:
3648 *INFORM = LUSOL_INFORM_NOPIVOT;
3649 if(LPRINT>=LUSOL_MSG_SINGULARITY)
3650 LUSOL_report(LUSOL, 0, "lu1fac error...\nTSP used but diagonal pivot could not be found\n");
3651
3652 /* Finalize and store output parameters. */
3653 x990:
3654 LUSOL->nelem = NELEM0;
3655 LUSOL->luparm[LUSOL_IP_SINGULARITIES] = NSING;
3656 LUSOL->luparm[LUSOL_IP_SINGULARINDEX] = JSING;
3657 LUSOL->luparm[LUSOL_IP_MINIMUMLENA] = MINLEN;
3658 LUSOL->luparm[LUSOL_IP_UPDATECOUNT] = 0;
3659 LUSOL->luparm[LUSOL_IP_RANK_U] = NRANK;
3660 LUSOL->luparm[LUSOL_IP_COLCOUNT_DENSE1] = NDENS1;
3661 LUSOL->luparm[LUSOL_IP_COLCOUNT_DENSE2] = NDENS2;
3662 LUSOL->luparm[LUSOL_IP_COLINDEX_DUMIN] = JUMIN;
3663 LUSOL->luparm[LUSOL_IP_COLCOUNT_L0] = NUML0;
3664 LUSOL->luparm[LUSOL_IP_ROWCOUNT_L0] = 0;
3665 LUSOL->luparm[LUSOL_IP_NONZEROS_L0] = LENL;
3666 LUSOL->luparm[LUSOL_IP_NONZEROS_U0] = LENU;
3667 LUSOL->luparm[LUSOL_IP_NONZEROS_L] = LENL;
3668 LUSOL->luparm[LUSOL_IP_NONZEROS_U] = LENU;
3669 LUSOL->luparm[LUSOL_IP_NONZEROS_ROW] = LROW;
3670 LUSOL->luparm[LUSOL_IP_MARKOWITZ_MERIT] = MERSUM;
3671 LUSOL->luparm[LUSOL_IP_TRIANGROWS_U] = NUTRI;
3672 LUSOL->luparm[LUSOL_IP_TRIANGROWS_L] = NLTRI;
3673 LUSOL->parmlu[LUSOL_RP_MAXELEM_A] = AMAX;
3674 LUSOL->parmlu[LUSOL_RP_MAXMULT_L] = LMAX;
3675 LUSOL->parmlu[LUSOL_RP_MAXELEM_U] = UMAX;
3676 LUSOL->parmlu[LUSOL_RP_MAXELEM_DIAGU] = DUMAX;
3677 LUSOL->parmlu[LUSOL_RP_MINELEM_DIAGU] = DUMIN;
3678 LUSOL->parmlu[LUSOL_RP_MAXELEM_TCP] = AKMAX;
3679 AGRWTH = AKMAX/(AMAX+LUSOL_SMALLNUM);
3680 UGRWTH = UMAX/(AMAX+LUSOL_SMALLNUM);
3681 if(TPP)
3682 GROWTH = UGRWTH;
3683 /* TRP or TCP or TSP */
3684 else
3685 GROWTH = AGRWTH;
3686 LUSOL->parmlu[LUSOL_RP_GROWTHRATE] = GROWTH;
3687
3688 LUSOL->luparm[LUSOL_IP_FTRANCOUNT] = 0;
3689 LUSOL->luparm[LUSOL_IP_BTRANCOUNT] = 0;
3690
3691 /* ------------------------------------------------------------------
3692 Set overall status variable.
3693 ------------------------------------------------------------------ */
3694 LUSOL->luparm[LUSOL_IP_INFORM] = *INFORM;
3695 if(*INFORM == LUSOL_INFORM_NOMEMLEFT)
3696 LUSOL_report(LUSOL, 0, "lu1fac error...\ninsufficient memory available\n");
3697
3698 /* ------------------------------------------------------------------
3699 Print statistics for the LU factors.
3700 ------------------------------------------------------------------ */
3701 NCP = LUSOL->luparm[LUSOL_IP_COMPRESSIONS_LU];
3702 CONDU = DUMAX/MAX(DUMIN,LUSOL_SMALLNUM);
3703 DINCR = (LENL+LENU)-LUSOL->nelem;
3704 DINCR = (DINCR*100)/MAX(DELEM,ONE);
3705 AVGMER = MERSUM;
3706 AVGMER = AVGMER/DM;
3707 NBUMP = LUSOL->m-NUTRI-NLTRI;
3708 if(LPRINT>=LUSOL_MSG_STATISTICS) {
3709 if(TPP) {
3710 LUSOL_report(LUSOL, 0, "Merit %g %d %d %d %g %d %d %g %g %d %d %d\n",
3711 AVGMER,LENL,LENL+LENU,NCP,DINCR,NUTRI,LENU,
3712 LTOL,UMAX,UGRWTH,NLTRI,NDENS1,LMAX);
3713 }
3714 else {
3715 LUSOL_report(LUSOL, 0, "Merit %s %g %d %d %d %g %d %d %g %g %d %d %d %g %g\n",
3716 LUSOL_pivotLabel(LUSOL),
3717 AVGMER,LENL,LENL+LENU,NCP,DINCR,NUTRI,LENU,
3718 LTOL,UMAX,UGRWTH,NLTRI,NDENS1,LMAX,AKMAX,AGRWTH);
3719 }
3720 LUSOL_report(LUSOL, 0, "bump%9d dense2%7d DUmax%g DUmin%g conDU%g\n",
3721 NBUMP,NDENS2,DUMAX,DUMIN,CONDU);
3722 }
3723 }
3724
3725
3726