1
2 /* Create a row-based version of L0.
3 This makes it possible to solve L0'x=h (btran) faster for sparse h,
4 since we only run down the columns of L0' (rows of LO) for which
5 the corresponding entry in h is non-zero. */
LU1L0(LUSOLrec * LUSOL,LUSOLmat ** mat,int * inform)6 MYBOOL LU1L0(LUSOLrec *LUSOL, LUSOLmat **mat, int *inform)
7 {
8 MYBOOL status = FALSE;
9 int K, L, LL, L1, L2, LENL0, NUML0, I;
10 int *lsumr;
11
12 /* Assume success */
13 *inform = LUSOL_INFORM_LUSUCCESS;
14
15 /* Check if there is anything worth doing */
16 if(mat == NULL)
17 return( status );
18 if(*mat != NULL)
19 LUSOL_matfree(mat);
20 NUML0 = LUSOL->luparm[LUSOL_IP_COLCOUNT_L0];
21 LENL0 = LUSOL->luparm[LUSOL_IP_NONZEROS_L0];
22 if((NUML0 == 0) || (LENL0 == 0) ||
23 (LUSOL->luparm[LUSOL_IP_ACCELERATION] == LUSOL_BASEORDER) ||
24 ((LUSOL->luparm[LUSOL_IP_ACCELERATION] & LUSOL_ACCELERATE_L0) == 0))
25 return( status );
26
27 /* Allocate temporary array */
28 lsumr = (int *) LUSOL_CALLOC((LUSOL->m+1), sizeof(*lsumr));
29 if(lsumr == NULL) {
30 *inform = LUSOL_INFORM_NOMEMLEFT;
31 return( status );
32 }
33
34 /* Compute non-zero counts by permuted row index (order is unimportant) */
35 K = 0;
36 L2 = LUSOL->lena;
37 L1 = L2-LENL0+1;
38 for(L = L1; L <= L2; L++) {
39 I = LUSOL->indc[L];
40 lsumr[I]++;
41 if(lsumr[I] == 1)
42 K++;
43 }
44 LUSOL->luparm[LUSOL_IP_ROWCOUNT_L0] = K;
45
46 /* Check if we should apply "smarts" before proceeding to the row matrix creation */
47 if((LUSOL->luparm[LUSOL_IP_ACCELERATION] & LUSOL_AUTOORDER) &&
48 ((REAL) LUSOL->luparm[LUSOL_IP_ROWCOUNT_L0] /
49 #if 0
50 LUSOL->luparm[LUSOL_IP_COLCOUNT_L0]
51 #else
52 LUSOL->m
53 #endif
54 > LUSOL->parmlu[LUSOL_RP_SMARTRATIO]))
55 goto Finish;
56
57 /* We are Ok to create the new matrix object */
58 *mat = LUSOL_matcreate(LUSOL->m, LENL0);
59 if(*mat == NULL) {
60 *inform = LUSOL_INFORM_NOMEMLEFT;
61 goto Finish;
62 }
63
64 /* Cumulate row counts to get vector offsets; first row is leftmost
65 (stick with Fortran array offset for consistency) */
66 (*mat)->lenx[0] = 1;
67 for(K = 1; K <= LUSOL->m; K++) {
68 (*mat)->lenx[K] = (*mat)->lenx[K-1] + lsumr[K];
69 lsumr[K] = (*mat)->lenx[K-1];
70 }
71
72 /* Map the matrix into row order by permuted index;
73 Note: The first permuted row is located leftmost in the array.
74 The column order is irrelevant, since the indeces will
75 refer to constant / resolved values of V[] during solve. */
76 L2 = LUSOL->lena;
77 L1 = L2-LENL0+1;
78 for(L = L1; L <= L2; L++) {
79 I = LUSOL->indc[L];
80 LL = lsumr[I]++;
81 (*mat)->a[LL] = LUSOL->a[L];
82 (*mat)->indr[LL] = LUSOL->indr[L];
83 (*mat)->indc[LL] = I;
84 }
85
86 /* Pack row starting positions, and set mapper from original index to packed */
87 I = 0;
88 for(L = 1; L <= LUSOL->m; L++) {
89 K = LUSOL->ip[L];
90 if((*mat)->lenx[K] > (*mat)->lenx[K-1]) {
91 I++;
92 (*mat)->indx[I] = K;
93 }
94 }
95
96 /* Confirm that everything went well */
97 status = TRUE;
98
99 /* Clean up */
100 Finish:
101 FREE(lsumr);
102 return( status );
103 }
104
105 /* Solve L0' v = v based on row-based version of L0, constructed by LU1L0 */
LU6L0T_v(LUSOLrec * LUSOL,LUSOLmat * mat,REAL V[],int NZidx[],int * INFORM)106 void LU6L0T_v(LUSOLrec *LUSOL, LUSOLmat *mat, REAL V[], int NZidx[], int *INFORM)
107 {
108 #ifdef DoTraceL0
109 REAL TEMP;
110 #endif
111 int LEN, K, KK, L, L1, NUML0;
112 REAL SMALL;
113 register REAL VPIV;
114 #if (defined LUSOLFastSolve) && !(defined DoTraceL0)
115 REAL *aptr;
116 int *jptr;
117 #else
118 int J;
119 #endif
120
121 NUML0 = LUSOL->luparm[LUSOL_IP_ROWCOUNT_L0];
122 SMALL = LUSOL->parmlu[LUSOL_RP_ZEROTOLERANCE];
123
124 /* Loop over the nz columns of L0' - from the end, going forward. */
125 for(K = NUML0; K > 0; K--) {
126 KK = mat->indx[K];
127 L = mat->lenx[KK];
128 L1 = mat->lenx[KK-1];
129 LEN = L - L1;
130 if(LEN == 0)
131 continue;
132 /* Get value of the corresponding active entry of V[] */
133 VPIV = V[KK];
134 /* Only process the column of L0' if the value of V[] is non-zero */
135 if(fabs(VPIV)>SMALL) {
136 /* ***** This loop could be coded specially. */
137 #if (defined LUSOLFastSolve) && !(defined DoTraceL0)
138 L--;
139 for(aptr = mat->a+L, jptr = mat->indr+L;
140 LEN > 0; LEN--, aptr--, jptr--)
141 V[*jptr] += VPIV * (*aptr);
142 #else
143 for(; LEN > 0; LEN--) {
144 L--;
145 J = mat->indr[L];
146 #ifndef DoTraceL0
147 V[J] += VPIV * mat->a[L];
148 #else
149 TEMP = V[J];
150 V[J] += VPIV * mat->a[L];
151 printf("V[%3d] = V[%3d] + L[%d,%d]*V[%3d]\n", J, J, KK,J, KK);
152 printf("%6g = %6g + %6g*%6g\n", V[J], TEMP, mat->a[L], VPIV);
153 #endif
154 }
155 #endif
156 }
157 #ifdef SetSmallToZero
158 else
159 V[KK] = 0;
160 #endif
161 }
162
163 }
164