1 /* bgmresr.c */
2
3 #include "../Iter.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 ---------------------------------------------------------------------
8 purpose -- to solve the system A X = B using the block GMRES
9 Arnoldi iteration method with right preconditioning
10
11 neqns -- # of equations
12 type -- must be SPOOLES_REAL or SPOOLES_COMPLEX
13 symmetryflag -- must be SPOOLES_SYMMETRIC or SPOOLES_NONSYMMETRIC
14 mtxA -- pointer to object for A
15 mtxM -- pointer to object for M, may be NULL
16 mtxX -- pointer to object for X, zero on input
17 mtxB -- pointer to object for B
18 maxnouter -- maximum number of outer iterations
19 maxninner -- maximum number of inner iterations
20 pninner -- the last # of inner iterations executed
21 pnouter -- the last # of outer iterations executed
22 convergetol -- tolerance for convergence
23 convergence achieved when ||r||_f <= convergetol * ||r_0||_f
24 msglvl -- message level
25 1 -- no output
26 2 -- iteration #, residual, ratio
27 3 -- scalar information
28 4 -- scalar and vector information
29 msgFile -- message file
30
31 return value ---
32 1 -- normal return
33 -1 -- neqns <= 0
34 -2 -- type invalid
35 -3 -- symmetryflag invalid
36 -4 -- mtxA NULL
37 -5 -- mtxX NULL
38 -6 -- mtxB NULL
39 -7 -- convergetol invalid
40 -8 -- maxnouter and/or maxninner invalid
41 -9 -- pnouter NULL
42 -10 -- pninner NULL
43 -11 -- msglvl > 0 and msgFile = NULL
44
45 created -- 98dec03, dkw
46 ---------------------------------------------------------------------
47 */
48 int
bgmresr(int neqns,int type,int symmetryflag,InpMtx * mtxA,FrontMtx * mtxM,DenseMtx * mtxX,DenseMtx * mtxB,int maxnouter,int maxninner,int * pnouter,int * pninner,double convergetol,int msglvl,FILE * msgFile)49 bgmresr (
50 int neqns,
51 int type,
52 int symmetryflag,
53 InpMtx *mtxA,
54 FrontMtx *mtxM,
55 DenseMtx *mtxX,
56 DenseMtx *mtxB,
57 int maxnouter,
58 int maxninner,
59 int *pnouter,
60 int *pninner,
61 double convergetol,
62 int msglvl,
63 FILE *msgFile
64 ) {
65 A2 *A2G, *A2Gi, *A2Gj, *A2H, *A2Hij, *A2Hkj, *A2H11, *A2H12,
66 *A2Vj, *A2Vk, *A2Z ;
67 DenseMtx *G, *Gj, *H, *Hij, *Hkj, *V, *Vi, *Vj, *Vk, *W, *Z;
68 double Hkk, Hik, initResNorm, ops, resnorm, t0, t1 ;
69 double cpus[9], minusone[2] = {-1.0, 0.0}, one[2] = {1.0, 0.0},
70 zero[2] = {0.0, 0.0} ;
71 double *col, *val, *Y ;
72 DV workDV ;
73 int converged, i, ii, iinner, jinner, jouter, kinner,
74 nrhs, nstep,
75 hijfrow, hijlrow, hijfcol, hijlcol, hkjfcol, hkjlcol,
76 iend, inc1, inc2, irow, istart, jcol, ncol, ndiag, nrow,
77 rc, vifcol, vilcol, vjfcol, vjlcol, vkfcol, vklcol ;
78 /*
79 ---------------
80 check the input
81 ---------------
82 */
83 MARKTIME(t0) ;
84 if ( neqns <= 0 ) {
85 fprintf(stderr, "\n error in bgmresr()"
86 "\n neqns = %d\n", neqns) ;
87 return(-1) ;
88 }
89 if ( type != SPOOLES_REAL ) {
90 fprintf(stderr, "\n error in bgmresr()"
91 "\n type = %d, must be SPOOLES_REAL\n", type) ;
92 return(-2) ;
93 }
94 if ( symmetryflag != SPOOLES_SYMMETRIC
95 && symmetryflag != SPOOLES_NONSYMMETRIC ) {
96 fprintf(stderr, "\n error in bgmresr()"
97 "\n symmetryflag = %d"
98 "\n must be SPOOLES_SYMMETRIC or SPOOLES_NONSYMMETRIC\n",
99 symmetryflag) ;
100 return(-3) ;
101 }
102 if ( mtxA == NULL ) {
103 fprintf(stderr, "\n error in bgmresr()"
104 "\n mtxA is NULL\n") ;
105 return(-4) ;
106 }
107 if ( mtxX == NULL ) {
108 fprintf(stderr, "\n error in bgmresr()"
109 "\n mtxX is NULL\n") ;
110 return(-5) ;
111 }
112 if ( mtxB == NULL ) {
113 fprintf(stderr, "\n error in bgmresr()"
114 "\n mtxB is NULL\n") ;
115 return(-6) ;
116 }
117 if ( convergetol < 0.0 ) {
118 fprintf(stderr, "\n error in bgmresr()"
119 "\n convergetol is %e\n", convergetol) ;
120 return(-7) ;
121 }
122 if ( maxnouter <= 0 || maxninner <= 0 ) {
123 fprintf(stderr, "\n error in bgmresr()"
124 "\n maxnouter %d, maxninner %d\n", maxnouter, maxninner) ;
125 return(-8) ;
126 }
127 if ( pnouter == NULL ) {
128 fprintf(stderr, "\n error in bgmresr()"
129 "\n pnouter is NULL\n") ;
130 return(-9) ;
131 }
132 if ( pninner == NULL ) {
133 fprintf(stderr, "\n error in bgmresr()"
134 "\n pninner is NULL\n") ;
135 return(-10) ;
136 }
137 if ( msglvl > 0 && msgFile == NULL ) {
138 fprintf(stderr, "\n error in bgmresr()"
139 "\n msglvl = %d and msgFile is NULL\n", msglvl) ;
140 return(-11) ;
141 }
142 if ( msglvl > 1 ) {
143 fprintf(msgFile,
144 "\n\n maximum # of outer iterations = %d"
145 "\n maximum # of inner iterations = %d",
146 maxnouter, maxninner) ;
147 fflush(msgFile) ;
148 }
149 /*
150 -------------------------------------------------
151 check for zero rhs (if B is zero, then X is zero)
152 -------------------------------------------------
153 */
154 initResNorm = DenseMtx_frobNorm(mtxB) ;
155 if ( msglvl > 1 ) {
156 fprintf(msgFile, "\n\n initial residual = %12.4e", initResNorm) ;
157 fflush(msgFile) ;
158 }
159 if ( initResNorm == 0.0 ) {
160 *pnouter = *pninner = 0 ;
161 return(1) ;
162 }
163 /*
164 ---------------------------
165 initialize the working data
166 ---------------------------
167 */
168 nrhs = mtxB->ncol ;
169 Y = DVinit(maxninner+1, 0.0) ;
170 DV_setDefaultFields(&workDV) ;
171
172 H = DenseMtx_new() ;
173 DenseMtx_init(H, SPOOLES_REAL, 0, 0,
174 (maxninner+1)*nrhs, maxninner*nrhs, 1, (maxninner+1)*nrhs) ;
175 DenseMtx_zero(H) ;
176
177 V = DenseMtx_new() ;
178 DenseMtx_init(V, SPOOLES_REAL, 0, 0, neqns, nrhs*(maxninner+1), 1, neqns) ;
179 DenseMtx_zero(V) ;
180
181 G = DenseMtx_new() ;
182 DenseMtx_init(G, SPOOLES_REAL, 0, 0,
183 (maxninner+1)*nrhs, nrhs, 1, (maxninner+1)*nrhs) ;
184 DenseMtx_zero(G) ;
185
186 Z = DenseMtx_new() ;
187 DenseMtx_init(Z, SPOOLES_REAL, 0, 0, neqns, nrhs, 1, neqns) ;
188 DenseMtx_zero(Z) ;
189
190 W = DenseMtx_new() ;
191 DenseMtx_init(W, SPOOLES_REAL, 0, 0, neqns, nrhs, 1, neqns) ;
192 DenseMtx_zero(W) ;
193
194 Vi = DenseMtx_new() ;
195 Vj = DenseMtx_new() ;
196 Vk = DenseMtx_new() ;
197 Gj = DenseMtx_new() ;
198 Hij = DenseMtx_new() ;
199 Hkj = DenseMtx_new() ;
200 A2G = A2_new() ;
201 A2H = A2_new() ;
202 A2Z = A2_new() ;
203 DenseMtx_setA2(Z, A2Z) ;
204 /*
205 ---------------------------
206 main loop of the iterations
207 ---------------------------
208 */
209 converged = 0 ;
210 for ( jouter = 0 ; jouter < maxnouter ; jouter++ ) {
211 A2Gi = A2_new() ;
212 A2Vj = A2_new() ;
213 DenseMtx_zero(H) ;
214 DenseMtx_zero(V) ;
215 DenseMtx_zero(G) ;
216 if ( msglvl > 2 ) {
217 fprintf(msgFile, "\n\n outer iteration %d", jouter) ;
218 fflush(msgFile) ;
219 }
220 /*
221 -----------------------------------------------------
222 compute the preconditioned residual and load into V_0
223 1. solve M Z = X
224 2. compute Z = A*Z
225 -----------------------------------------------------
226 */
227 if ( msglvl > 3 ) {
228 fprintf(msgFile, "\n\n B") ;
229 DenseMtx_writeForHumanEye(mtxB, msgFile) ;
230 fflush(msgFile) ;
231 }
232 if ( jouter == 0 ) {
233 for ( i = 0 ; i < nrhs ; i++ )
234 DenseMtx_colCopy(Z, i, mtxB, i) ;
235 }else{
236 if ( mtxM != NULL ) {
237 FrontMtx_solve(mtxM, W, mtxX, mtxM->manager,
238 cpus, 0, msgFile) ;
239 }
240 for ( i = 0 ; i < nrhs ; i++ )
241 DenseMtx_colCopy(Z, i, mtxB, i) ;
242 if ( symmetryflag == SPOOLES_SYMMETRIC ) {
243 InpMtx_sym_mmm(mtxA, Z, minusone, W) ;
244 } else {
245 InpMtx_nonsym_mmm(mtxA, Z, minusone, W) ;
246 }
247 }
248 /*
249 ---------------------------------------------------
250 Compute the QR factorization of the initial
251 residual to get our starting orthogonal vectors V_0
252 and our top block in G
253 ---------------------------------------------------
254 */
255 ops = A2_QRreduce(A2Z, &workDV, msglvl, msgFile) ;
256
257 rc = DenseMtx_initAsSubmatrix(Vj, V, 0, neqns-1, 0, nrhs-1) ;
258 DenseMtx_setA2(Vj, A2Vj) ;
259 A2_computeQ(A2Vj, A2Z, &workDV, msglvl, msgFile) ;
260
261 DenseMtx_initAsSubmatrix(Gj, G, 0, 2*nrhs-1, 0, nrhs-1) ;
262 DenseMtx_setA2(Gj, A2Gi) ;
263
264 nrow = A2Z->n1 ;
265 ncol = A2Z->n2 ;
266 inc1 = A2Z->inc1 ;
267 inc2 = A2Z->inc2 ;
268 if ( nrow >= ncol ) {
269 ndiag = ncol ;
270 } else {
271 ndiag = nrow ;
272 }
273 for ( jcol = 0, col = A2Z->entries, val = A2Gi->entries ;
274 jcol < ncol ;
275 jcol++, col += inc2, val += inc2 ) {
276 istart = 0 ;
277 iend = (jcol < ndiag) ? jcol : ndiag - 1 ;
278 for ( irow = istart, ii = irow*inc1 ;
279 irow <= iend ;
280 irow++, ii += inc1 ) {
281 val[ii] = col[ii] ;
282 }
283 }
284 resnorm = DenseMtx_frobNorm(G) ;
285 if ( msglvl > 2 ) {
286 fprintf(msgFile, "\n\n after QRreduce A2Vj =") ;
287 fprintf(msgFile, "\n\n V_%d = ", jouter) ;
288 A2_writeForHumanEye(A2Vj, msgFile) ;
289 fprintf(msgFile, "\n\n G_%d%d = ", jouter, jouter) ;
290 A2_writeForHumanEye(A2Gi, msgFile) ;
291 fflush(msgFile) ;
292 }
293 if ( msglvl > 2 ) {
294 fprintf(msgFile, "\n\n ||G||_f = %12.4e", resnorm) ;
295 fflush(msgFile) ;
296 }
297 if ( resnorm == 0.0 ) {
298 break ;
299 }
300 if ( jouter == 0 ) {
301 initResNorm = resnorm ;
302 }
303 /*
304 ------------------------------------------------------------
305 loop over the inner gmres iterations using Arnoldi iteration
306 ------------------------------------------------------------
307 */
308 nstep = maxninner - 1 ;
309 for ( jinner = 0 ; jinner < maxninner ; jinner++ ) {
310 A2Gj = A2_new() ;
311 A2Hij = A2_new() ;
312 A2Hkj = A2_new() ;
313 if ( msglvl > 2 ) {
314 fprintf(msgFile, "\n\n inner iteration %d", jinner) ;
315 fflush(msgFile) ;
316 }
317 /*
318 ----------------------------
319 compute W = A * M^{-1} * V_j
320 ----------------------------
321 */
322 DenseMtx_zero(Z) ;
323 vjfcol = jinner*nrhs ;
324 vjlcol = vjfcol+nrhs-1 ;
325 DenseMtx_initAsSubmatrix(Vj, V, 0, neqns-1, vjfcol, vjlcol) ;
326 if ( mtxM != NULL ) {
327 FrontMtx_solve(mtxM, W, Vj, mtxM->manager,
328 cpus, 0, msgFile) ;
329 }
330 if ( msglvl > 2 ) {
331 fprintf(msgFile, "\n\n ||M^{-1}*V_%d||_f = %12.4e",
332 jinner, DenseMtx_frobNorm(W)) ;
333 fflush(msgFile) ;
334 }
335 if ( symmetryflag == SPOOLES_SYMMETRIC ) {
336 InpMtx_sym_mmm(mtxA, Z, one, W) ;
337 } else {
338 InpMtx_nonsym_mmm(mtxA, Z, one, W) ;
339 }
340 if ( msglvl > 2 ) {
341 fprintf(msgFile, "\n\n ||A*M^{-1}V_%d||_f = %12.4e",
342 jinner, DenseMtx_frobNorm(Z)) ;
343 fflush(msgFile) ;
344 }
345 if ( msglvl > 3 ) {
346 fprintf(msgFile, "\n\n A*M^{-1}*V_%d", jinner) ;
347 DenseMtx_writeForHumanEye(Z, msgFile) ;
348 fflush(msgFile) ;
349 }
350
351 for ( iinner = 0 ; iinner <= jinner ; iinner++ ) {
352 vifcol = iinner*nrhs ;
353 vilcol = vifcol+nrhs-1 ;
354 DenseMtx_initAsSubmatrix(Vi, V, 0, neqns-1, vifcol, vilcol) ;
355 hijfrow = vifcol ;
356 hijlrow = vilcol ;
357 hijfcol = jinner*nrhs ;
358 hijlcol = hijfcol+nrhs-1 ;
359 DenseMtx_initAsSubmatrix(Hij, H, hijfrow, hijlrow, hijfcol, hijlcol) ;
360 DenseMtx_mmm("t", "n", zero, Hij, one, Vi, Z) ;
361 DenseMtx_mmm("n", "n", one, Z, minusone, Vi, Hij) ;
362 if ( msglvl > 2 ) {
363 fprintf(msgFile, "\n\n H_%d%d = ", iinner, jinner) ;
364 DenseMtx_writeForHumanEye(Hij, msgFile) ;
365 fflush(msgFile) ;
366 }
367 }
368 if ( msglvl > 3 ) {
369 fprintf(msgFile, "\n\n Wj after Arnoldi iteration") ;
370 DenseMtx_writeForHumanEye(Z, msgFile) ;
371 fflush(msgFile) ;
372 }
373 /*
374 -----------------------------------------------
375 compute V_jinner+1 and H_jinner+1,jinner blocks
376 using QR decomposition of W_j then
377 copy R from A2Z(1:ncolA,1:ncolA) into H_jinner+1,jinner
378 -----------------------------------------------
379 */
380 ops = A2_QRreduce(A2Z, &workDV, msglvl, msgFile) ;
381 vkfcol = vjlcol+1 ;
382 vklcol = vkfcol+nrhs-1 ;
383 DenseMtx_initAsSubmatrix(Vk, V, 0, neqns-1, vkfcol, vklcol) ;
384 A2Vk = A2_new() ;
385 DenseMtx_setA2(Vk, A2Vk) ;
386 A2_computeQ(A2Vk, A2Z, &workDV, msglvl, msgFile) ;
387 A2_free(A2Vk) ;
388
389 hijfrow = hijlrow+1 ;
390 hijlrow = hijfrow+nrhs-1 ;
391 DenseMtx_initAsSubmatrix(Hkj, H, hijfrow, hijlrow, hijfcol, hijlcol) ;
392 DenseMtx_setA2(Hkj, A2Hkj) ;
393 nrow = A2Z->n1 ;
394 ncol = A2Z->n2 ;
395 inc1 = A2Z->inc1 ;
396 inc2 = A2Z->inc2 ;
397 if ( nrow >= ncol ) {
398 ndiag = ncol ;
399 } else {
400 ndiag = nrow ;
401 }
402 for ( jcol = 0, col = A2Z->entries, val = A2Hkj->entries ;
403 jcol < ncol ;
404 jcol++, col += inc2, val += inc2 ) {
405 istart = 0 ;
406 iend = (jcol < ndiag) ? jcol : ndiag - 1 ;
407 for ( irow = istart, ii = irow*inc1 ;
408 irow <= iend ;
409 irow++, ii += inc1 ) {
410 val[ii] = col[ii] ;
411 }
412 }
413 if ( msglvl > 2 ) {
414 fprintf(msgFile, "\n\n V_%d = ", jinner+1) ;
415 DenseMtx_writeForHumanEye(Vk, msgFile) ;
416 fprintf(msgFile, "\n\n H_%d%d = ", jinner+1, jinner) ;
417 DenseMtx_writeForHumanEye(Hij, msgFile) ;
418 fflush(msgFile) ;
419 }
420 /*
421 -------------------------------------------------
422 Before starting the next iteration, triangulate
423 the Hessenberg matrix. To do this we must first
424 apply the previous triangulating transformations
425 to the newly constructed columns and then
426 triangulate the new portion of the Hessenberg
427 matrix. When we compute these new transformations
428 we need to also apply them to the corresponding
429 portion of G to keep track of the norm of the
430 current residual.
431
432 First apply the previous transformations
433 -------------------------------------------------
434 */
435 hijfrow = 0 ;
436 hijlrow = 2*nrhs-1 ;
437 hkjfcol = 0 ;
438 hkjlcol = nrhs-1 ;
439 for ( iinner = 0 ; iinner <= jinner-1 ; iinner++ ) {
440 A2H11 = A2_new() ;
441 A2H12 = A2_new() ;
442 DenseMtx_initAsSubmatrix(Hij, H, hijfrow, hijlrow, hkjfcol, hkjlcol) ;
443 DenseMtx_setA2(Hij, A2H11) ;
444 DenseMtx_initAsSubmatrix(Hkj, H, hijfrow, hijlrow, hijfcol, hijlcol) ;
445 DenseMtx_setA2(Hkj, A2H12) ;
446 A2_applyQT(A2H12, A2H11, A2H12, &workDV, msglvl, msgFile) ;
447 hijfrow += nrhs ;
448 hijlrow += nrhs ;
449 hkjfcol += nrhs ;
450 hkjlcol += nrhs ;
451 A2_free(A2H11) ;
452 A2_free(A2H12) ;
453 }
454 /*
455 -------------------------------------------------
456 Compute the QR factorization of the diagonal and
457 subdiagonal blocks of the block Hessenberg matrix
458 -------------------------------------------------
459 */
460 DenseMtx_initAsSubmatrix(Hij, H, hijfrow, hijlrow, hijfcol, hijlcol) ;
461 DenseMtx_setA2(Hij, A2Hij) ;
462 ops = A2_QRreduce(A2Hij, &workDV, msglvl, msgFile) ;
463 /*
464 ----------------------------------------------------------
465 apply the transformation to the corresponding section of G
466 ----------------------------------------------------------
467 */
468 DenseMtx_initAsSubmatrix(Gj, G, hijfrow, hijlrow, 0, nrhs-1) ;
469 DenseMtx_setA2(Gj, A2Gj) ;
470 A2_applyQT(A2Gj, A2Hij, A2Gj, &workDV, msglvl, msgFile) ;
471 DenseMtx_initAsSubmatrix(Gj, G, hijfrow+nrhs, hijlrow, 0, nrhs-1) ;
472 resnorm = DenseMtx_frobNorm(Gj) ;
473 if ( msglvl > 1 ) {
474 fprintf(msgFile, "\n outer %d, inner %d, res = %24.16e",
475 jouter, jinner, resnorm) ;
476 fflush(msgFile) ;
477 }
478 /*
479 -----------------
480 check convergence
481 -----------------
482 */
483 if ( resnorm <= convergetol * initResNorm ) {
484 if ( msglvl > 1 ) {
485 fprintf(msgFile, "\n BGMRESR convergence has been achieved") ;
486 fprintf(msgFile, "\n ***convergetol = %24.16e", convergetol) ;
487 fprintf(msgFile, "\n ***initResNorm = %24.16e", initResNorm) ;
488 fflush(msgFile) ;
489 }
490 /*
491 ------------------------------------------------
492 convergence has been achieved, break out of loop
493 ------------------------------------------------
494 */
495 nstep = jinner ;
496 converged = 1 ;
497 A2_free(A2Gj) ;
498 A2_free(A2Hij) ;
499 A2_free(A2Hkj) ;
500 break ;
501 }
502 A2_free(A2Gj) ;
503 A2_free(A2Hij) ;
504 A2_free(A2Hkj) ;
505 }
506 /*
507 ----------------------------------
508 for each rhs copy Gi into vector Y
509 solve Y := H^{-1} Y
510 and store Y back into Gi
511 ----------------------------------
512 */
513 DenseMtx_setA2(H, A2H) ;
514 DenseMtx_setA2(G, A2G) ;
515 if ( msglvl > 2 ) {
516 fprintf(msgFile, "\n\n H") ;
517 A2H->n1 = nstep+1 ;
518 A2H->n2 = nstep+1 ;
519 A2_writeForHumanEye(A2H, msgFile) ;
520 A2G->n1 = nstep+1 ;
521 fprintf(msgFile, "\n\n G") ;
522 A2_writeForHumanEye(A2G, msgFile) ;
523 A2H->n1 = maxninner+1 ;
524 A2H->n2 = maxninner+1 ;
525 fprintf(msgFile, "\n\n before solve, Y") ;
526 fprintf(msgFile, "\n\n before solve ykinner = %24.16e ",Y[7]) ;
527 DVfprintf(msgFile, nstep+1, Y) ;
528 fflush(msgFile) ;
529 }
530 for ( i = 0 ; i < nrhs ; i++ ) {
531 A2_extractColumn(A2G, Y, i) ;
532 for ( kinner = nstep ; kinner >= 0 ; kinner-- ) {
533 A2_realEntry(A2H, kinner, kinner, &Hkk) ;
534 Y[kinner] = Y[kinner] / Hkk ;
535 for ( iinner = 0 ; iinner < kinner ; iinner++ ) {
536 A2_realEntry(A2H, iinner, kinner, &Hik) ;
537 Y[iinner] -= Y[kinner] * Hik ;
538 }
539 }
540 if ( msglvl > 2 ) {
541 fprintf(msgFile, "\n\n after solve, Y") ;
542 DVfprintf(msgFile, nstep+1, Y) ;
543 fflush(msgFile) ;
544 }
545 A2_setColumn(A2G, Y, i) ;
546 }
547 /*
548 -------------------------------------
549 compute the new solution X := X + V*Y
550 -------------------------------------
551 */
552 DenseMtx_mmm("n", "n", one, mtxX, one, V, G) ;
553 if ( msglvl > 2 ) {
554 resnorm = DenseMtx_frobNorm(mtxX) ;
555 fprintf(msgFile, "\n ||X||_f = %24.16e", resnorm) ;
556 fflush(msgFile) ;
557 }
558 if ( msglvl > 3 ) {
559 fprintf(msgFile, "\n\n X") ;
560 DenseMtx_writeForHumanEye(mtxX, msgFile) ;
561 }
562 if ( converged == 1 ) {
563 break ;
564 }
565 A2_free(A2Gi) ;
566 A2_free(A2Vj) ;
567 }
568 /*
569 ------------------
570 final solve MX = Y
571 ------------------
572 */
573 if ( mtxM != NULL ) {
574 FrontMtx_solve(mtxM, mtxX, mtxX, mtxM->manager,
575 cpus, 0, msgFile) ;
576 }
577 if ( msglvl > 2 ) {
578 resnorm = DenseMtx_frobNorm(mtxX) ;
579 fprintf(msgFile, "\n ***after solve ||X||_f = %24.16e", resnorm) ;
580 }
581 /*
582 -----------
583 return info
584 -----------
585 */
586 *pnouter = jouter + 1 ;
587 *pninner = jinner + 1 ;
588 if ( msglvl > 2 ) {
589 fprintf(msgFile, "\n\n before leaving, *pnouter = %d, *pinner = %d",
590 *pnouter, *pninner) ;
591 fflush(msgFile) ;
592 }
593 /*
594 ------------------------
595 free the working storage
596 ------------------------
597 */
598 A2_free(A2G) ;
599 A2_free(A2H) ;
600 A2_free(A2Z) ;
601 A2_free(A2Gi) ;
602 A2_free(A2Vj) ;
603 DenseMtx_free(Gj) ;
604 DenseMtx_free(G) ;
605 DenseMtx_free(Hij) ;
606 DenseMtx_free(Hkj) ;
607 DenseMtx_free(H) ;
608 DenseMtx_free(Vi) ;
609 DenseMtx_free(Vj) ;
610 DenseMtx_free(Vk) ;
611 DenseMtx_free(V) ;
612 DenseMtx_free(Z) ;
613 DenseMtx_free(W) ;
614 DV_clearData(&workDV) ;
615 DVfree(Y) ;
616
617 return(1) ; }
618
619 /*--------------------------------------------------------------------*/
620