1 /* CalculiX - A 3-dimensional finite element program */
2 /* Copyright (C) 1998-2021 Guido Dhondt */
3
4 /* This program is free software; you can redistribute it and/or */
5 /* modify it under the terms of the GNU General Public License as */
6 /* published by the Free Software Foundation(version 2); */
7 /* */
8
9 /* This program is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU General Public License for more details. */
13
14 /* You should have received a copy of the GNU General Public License */
15 /* along with this program; if not, write to the Free Software */
16 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
17
18 #include <stdio.h>
19 #include <math.h>
20 #include <stdlib.h>
21 #include <string.h>
22
23 #ifdef SPOOLES
24 #include <misc.h>
25 #include <FrontMtx.h>
26 #include <SymbFac.h>
27 #endif
28
29 #include "CalculiX.h"
30
31 #define min(a,b) ((a) <= (b) ? (a) : (b))
32 #define max(a,b) ((a) >= (b) ? (a) : (b))
33
cascadefem(ITG * ipompc,double ** coefmpcp,ITG ** nodempcp,ITG * nmpc,ITG * mpcfree,ITG * nodeboun,ITG * ndirboun,ITG * nboun,ITG * ikmpc,ITG * ilmpc,ITG * ikboun,ITG * ilboun,ITG * mpcend,ITG * mpcmult,char * labmpc,ITG * nk,ITG * memmpc_,ITG * icascade,ITG * maxlenmpc,ITG * callfrommain,ITG * iperturb,ITG * ithermal)34 void cascadefem(ITG *ipompc, double **coefmpcp, ITG **nodempcp, ITG *nmpc,
35 ITG *mpcfree, ITG *nodeboun, ITG *ndirboun, ITG*nboun, ITG*ikmpc,
36 ITG *ilmpc, ITG *ikboun, ITG *ilboun, ITG *mpcend, ITG *mpcmult,
37 char *labmpc, ITG *nk, ITG *memmpc_, ITG *icascade, ITG *maxlenmpc,
38 ITG *callfrommain, ITG *iperturb, ITG *ithermal){
39
40 /* detects cascaded mpc's and decascades them; checks multiple
41 occurrence of the same dependent DOF's in different mpc/spc's
42
43 data structure of ipompc,coefmpc,nodempc:
44 for each mpc, e.g. i,
45 -the nodes are stored in nodempc(1,ipompc(i)),
46 nodempc(1,nodempc(3,ipompc(i))),
47 nodempc(1,nodempc(3,nodempc(3,ipompc(i))))... till
48 nodempc(3,nodempc(3,nodempc(3,.......))))))=0;
49 -the corresponding directions in nodempc(2,ipompc(i)),
50 nodempc(2,nodempc(3,ipompc(i))),.....
51 -the corresponding coefficient in coefmpc(ipompc(i)),
52 coefmpc(nodempc(3,ipompc(i))),.....
53 the mpc is written as a(1)u(i1,j1)+a(2)u(i2,j2)+...
54 +....a(k)u(ik,jk)=0, the first term is the dependent term,
55 the others are independent, at least after execution of the
56 present routine. The mpc's must be homogeneous, otherwise a
57 error message is generated and the program stops. */
58
59 ITG i,j,index,id,idof,nterm,idepend,*nodempc=NULL,
60 ispooles,iexpand,ichange,indexold,
61 mpc,indexnew,index1,index2,index1old,index2old,*jmpc=NULL,nl;
62
63 double coef,*coefmpc=NULL;
64
65 #ifdef SPOOLES
66
67 ITG irow,icolumn,node,idir,irownl,icolnl,*ipointer=NULL,*icoef=NULL,
68 ifree,*indepdof=NULL,nindep;
69
70 double *xcoef=NULL,b;
71
72 DenseMtx *mtxB, *mtxX ;
73 Chv *rootchv ;
74 ChvManager *chvmanager ;
75 SubMtxManager *mtxmanager ;
76 FrontMtx *frontmtx ;
77 InpMtx *mtxA ;
78 double tau = 100.;
79 double cpus[10] ;
80 ETree *frontETree ;
81 FILE *msgFile ;
82 Graph *graph ;
83 ITG jrhs, msglvl=0, nedges,error,
84 nent, neqns, nrhs, pivotingflag=1, seed=389,
85 symmetryflag=2, type=1,maxdomainsize,maxzeros,maxsize;
86 ITG *oldToNew ;
87 ITG stats[20] ;
88 IV *newToOldIV, *oldToNewIV ;
89 IVL *adjIVL, *symbfacIVL ;
90 #endif
91
92 nodempc=*nodempcp;
93 coefmpc=*coefmpcp;
94
95 /* for(i=0;i<*nmpc;i++){
96 j=i+1;
97 FORTRAN(writempc,(ipompc,nodempc,coefmpc,labmpc,&j));
98 }*/
99
100 NNEW(jmpc,ITG,*nmpc);
101 idepend=0;
102
103 /* check whether a node is used as a dependent node in a MPC
104 and in a SPC */
105
106 for(i=0;i<*nmpc;i++){
107 if(*nboun>0){
108 FORTRAN(nident,(ikboun,&ikmpc[i],nboun,&id));}
109 else{id=0;}
110 if(id>0){
111 if(ikboun[id-1]==ikmpc[i]){
112 printf("*ERROR in cascade: the DOF corresponding to \n node %d in direction %d is detected on the \n dependent side of a MPC and a SPC\n",
113 (ikmpc[i])/8+1,ikmpc[i]-8*((ikmpc[i])/8));
114 FORTRAN(stop,());
115 }
116 }
117 }
118
119 /* check whether there are user mpc's: in user MPC's the
120 dependent DOF can change, however, the number of terms
121 cannot change */
122
123 for(i=0;i<*nmpc;i++){
124
125 /* linear mpc */
126
127 /* because of the next line the size of field labmpc
128 has to be defined as 20*nmpc+1: without "+1" an
129 undefined field is accessed */
130
131 if((strcmp1(&labmpc[20*i]," ")==0) ||
132 (strcmp1(&labmpc[20*i],"CYCLIC")==0) ||
133 (strcmp1(&labmpc[20*i],"SUBCYCLIC")==0)||
134 (strcmp1(&labmpc[20*i],"CONTACT")==0)||
135 (*iperturb<2)) jmpc[i]=0;
136
137 /* nonlinear mpc */
138
139 else if((strcmp1(&labmpc[20*i],"RIGID")==0) ||
140 (strcmp1(&labmpc[20*i],"KNOT")==0) ||
141 (strcmp1(&labmpc[20*i],"PLANE")==0) ||
142 (strcmp1(&labmpc[20*i],"STRAIGHT")==0)||
143 (strcmp1(&labmpc[20*i],"ISOCHORIC")==0)) jmpc[i]=1;
144
145 /* user mpc */
146
147 else{
148 jmpc[i]=1;
149 if(*icascade==0) *icascade=1;
150 }
151 }
152
153 /* decascading */
154
155 ispooles=0;
156
157 /* decascading using simple substitution */
158
159 do{
160 ichange=0;
161 for(i=0;i<*nmpc;i++){
162 if(jmpc[i]==1) nl=1;
163 else nl=0;
164 iexpand=0;
165 index=nodempc[3*ipompc[i]-1];
166 if(index==0) continue;
167 do{
168 idof=(nodempc[3*index-3]-1)*8+nodempc[3*index-2];
169 FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
170 if((id>0)&&(ikmpc[id-1]==idof)){
171
172 /* a term on the independent side of the MPC is
173 detected as dependent node in another MPC */
174
175 indexold=nodempc[3*index-1];
176 coef=coefmpc[index-1];
177 mpc=ilmpc[id-1];
178
179 /* no expansion if there is a dependence of a
180 nonlinear MPC on another linear or nonlinear MPC
181 and the call is from main */
182
183 if((jmpc[mpc-1]==1)||(nl==1)){
184 *icascade=2;
185 if(idepend==0){
186 printf("*INFO in cascade: linear MPCs and\n");
187 printf(" nonlinear MPCs depend on each other\n");
188 printf(" common node: %d in direction %d\n\n",nodempc[3*index-3],nodempc[3*index-2]);
189 idepend=1;}
190 if(*callfrommain==1){
191 index=nodempc[3*index-1];
192 if(index!=0) continue;
193 else break;}
194 }
195
196 /* printf("*INFO in cascade: DOF %d of node %d is expanded\n",
197 nodempc[3*index-2],nodempc[3*index-3]);*/
198
199 /* collecting terms corresponding to the same DOF */
200
201 index1=ipompc[i];
202 do{
203 index2old=index1;
204 index2=nodempc[3*index1-1];
205 if(index2==0) break;
206 do{
207 if((nodempc[3*index1-3]==nodempc[3*index2-3])&&
208 (nodempc[3*index1-2]==nodempc[3*index2-2])){
209 coefmpc[index1-1]+=coefmpc[index2-1];
210 nodempc[3*index2old-1]=nodempc[3*index2-1];
211 nodempc[3*index2-1]=*mpcfree;
212 *mpcfree=index2;
213 index2=nodempc[3*index2old-1];
214 if(index2==0) break;
215 }
216 else{
217 index2old=index2;
218 index2=nodempc[3*index2-1];
219 if(index2==0) break;
220 }
221 }while(1);
222 index1=nodempc[3*index1-1];
223 if(index1==0) break;
224 }while(1);
225
226 /* check for zero coefficients on the dependent side */
227
228 index1=ipompc[i];
229 if(fabs(coefmpc[index1-1])<1.e-10){
230 printf("*ERROR in cascade: zero coefficient on the\n");
231 printf(" dependent side of an equation\n");
232 printf(" dependent node: %d",nodempc[3*index1-3]);
233 FORTRAN(stop,());
234 }
235
236 ichange=1;iexpand=1;
237 if((strcmp1(&labmpc[20*i]," ")==0)&&
238 (strcmp1(&labmpc[20*(mpc-1)],"CYCLIC")==0))
239 strcpy1(&labmpc[20*i],"SUBCYCLIC",9);
240 indexnew=ipompc[mpc-1];
241 coef=-coef/coefmpc[indexnew-1];
242 indexnew=nodempc[3*indexnew-1];
243 do{
244 coefmpc[index-1]=coef*coefmpc[indexnew-1];
245 nodempc[3*index-3]=nodempc[3*indexnew-3];
246 nodempc[3*index-2]=nodempc[3*indexnew-2];
247 indexnew=nodempc[3*indexnew-1];
248 if(indexnew!=0){
249 nodempc[3*index-1]=*mpcfree;
250 index=*mpcfree;
251 *mpcfree=nodempc[3**mpcfree-1];
252 if(*mpcfree==0){
253 *mpcfree=*memmpc_+1;
254 nodempc[3*index-1]=*mpcfree;
255 *memmpc_=(ITG)(1.1**memmpc_);
256 printf("*INFO in cascade: reallocating nodempc; new size = %d\n\n",*memmpc_);
257 RENEW(nodempc,ITG,3**memmpc_);
258 RENEW(coefmpc,double,*memmpc_);
259 for(j=*mpcfree;j<*memmpc_;j++){
260 nodempc[3*j-1]=j+1;
261 }
262 nodempc[3**memmpc_-1]=0;
263 }
264 continue;
265 }
266 else{
267 nodempc[3*index-1]=indexold;
268 break;
269 }
270 }while(1);
271 break;
272 }
273 else{
274 index=nodempc[3*index-1];
275 if(index!=0) continue;
276 else break;
277 }
278 }while(1);
279 if(iexpand==0) continue;
280
281 /* one term of the mpc was expanded
282 collecting terms corresponding to the same DOF */
283
284 index1=ipompc[i];
285 do{
286 index2old=index1;
287 index2=nodempc[3*index1-1];
288 if(index2==0) break;
289 do{
290 if((nodempc[3*index1-3]==nodempc[3*index2-3])&&
291 (nodempc[3*index1-2]==nodempc[3*index2-2])){
292 coefmpc[index1-1]+=coefmpc[index2-1];
293 nodempc[3*index2old-1]=nodempc[3*index2-1];
294 nodempc[3*index2-1]=*mpcfree;
295 *mpcfree=index2;
296 index2=nodempc[3*index2old-1];
297 if(index2==0) break;
298 }
299 else{
300 index2old=index2;
301 index2=nodempc[3*index2-1];
302 if(index2==0) break;
303 }
304 }while(1);
305 index1=nodempc[3*index1-1];
306 if(index1==0) break;
307 }while(1);
308
309 /* check for zero coefficients on the dependent and
310 independent side */
311
312 index1=ipompc[i];
313 index1old=0;
314 do {
315 if(fabs(coefmpc[index1-1])<1.e-10){
316 if(index1old==0){
317 printf("*ERROR in cascade: zero coefficient on the\n");
318 printf(" dependent side of an equation\n");
319 printf(" dependent node: %d",nodempc[3*index1-3]);
320 FORTRAN(stop,());
321 }
322 else{
323 nodempc[3*index1old-1]=nodempc[3*index1-1];
324 nodempc[3*index1-1]=*mpcfree;
325 *mpcfree=index1;
326 index1=nodempc[3*index1old-1];
327 }
328 }
329 else{
330 index1old=index1;
331 index1=nodempc[3*index1-1];
332 }
333 if(index1==0) break;
334 }while(1);
335 }
336 if(ichange==0) break;
337 }while(1);
338
339 /* decascading using spooles */
340
341 #ifdef SPOOLES
342 if((*icascade==1)&&(ispooles==1)){
343 if ( (msgFile = fopen("spooles.out", "a")) == NULL ) {
344 fprintf(stderr, "\n fatal error in spooles.c"
345 "\n unable to open file spooles.out\n") ;
346 }
347 NNEW(ipointer,ITG,7**nk);
348 NNEW(indepdof,ITG,7**nk);
349 NNEW(icoef,ITG,2**memmpc_);
350 NNEW(xcoef,double,*memmpc_);
351 ifree=0;
352 nindep=0;
353
354 for(i=*nmpc-1;i>-1;i--){
355 index=ipompc[i];
356 while(1){
357 idof=8*(nodempc[3*index-3]-1)+nodempc[3*index-2]-1;
358
359 /* check whether idof is a independent dof which has not yet been
360 stored in indepdof */
361
362 FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
363 if((id==0)||(ikmpc[id-1]!=idof)){
364 FORTRAN(nident,(indepdof,&idof,&nindep,&id));
365 if((id==0)||(indepdof[id-1]!=idof)){
366 for(j=nindep;j>id;j--){
367 indepdof[j]=indepdof[j-1];
368 }
369 indepdof[id]=idof;
370 nindep++;
371 }
372 }
373
374 icoef[2*ifree]=i+1;
375 icoef[2*ifree+1]=ipointer[idof];
376 xcoef[ifree]=coefmpc[index-1];
377 ipointer[idof]=++ifree;
378 index=nodempc[3*index-1];
379 if(index==0) break;
380 }
381 }
382
383 /* filling the left hand side */
384
385 nent=*memmpc_;
386 neqns=*nmpc;
387 mtxA = InpMtx_new() ;
388 InpMtx_init(mtxA, INPMTX_BY_ROWS, type, nent, neqns) ;
389
390 for(i=0;i<*nmpc;i++){
391 idof=ikmpc[i];
392 icolumn=ilmpc[i]-1;
393 if(strcmp1(&labmpc[20*icolumn],"RIGID")==0) icolnl=1;
394 else icolnl=0;
395 index=ipointer[idof-1];
396 while(1){
397 irow=icoef[2*index-2]-1;
398 if(irow!=icolumn){
399 if(strcmp1(&labmpc[20*irow],"RIGID")==0)irownl=1;
400 else irownl=0;
401 if((irownl==1)||(icolnl==1)){
402 *icascade=2;
403 InpMtx_free(mtxA);
404 printf("*ERROR in cascade: linear and nonlinear MPCs depend on each other");
405 FORTRAN(stop,());
406 }
407 }
408 if((strcmp1(&labmpc[20*irow]," ")==0)&&
409 (strcmp1(&labmpc[20*icolumn],"CYCLIC")==0)){
410 strcpy1(&labmpc[20*irow],"SUBCYCLIC",9);}
411 coef=xcoef[index-1];
412 InpMtx_inputRealEntry(mtxA,irow,icolumn,coef);
413 index=icoef[2*index-1];
414 if(index==0) break;
415 }
416 ipointer[idof-1]=0;
417 }
418
419 InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ;
420 if ( msglvl > 1 ) {
421 fprintf(msgFile, "\n\n input matrix") ;
422 InpMtx_writeForHumanEye(mtxA, msgFile) ;
423 fflush(msgFile) ;
424 }
425 /*--------------------------------------------------------------------*/
426 /*
427 -------------------------------------------------
428 STEP 2 : find a low-fill ordering
429 (1) create the Graph object
430 (2) order the graph using multiple minimum degree
431 -------------------------------------------------
432 */
433 graph = Graph_new() ;
434 adjIVL = InpMtx_fullAdjacency(mtxA) ;
435 nedges = IVL_tsize(adjIVL) ;
436 Graph_init2(graph, 0, neqns, 0, nedges, neqns, nedges, adjIVL,
437 NULL, NULL) ;
438 if ( msglvl > 1 ) {
439 fprintf(msgFile, "\n\n graph of the input matrix") ;
440 Graph_writeForHumanEye(graph, msgFile) ;
441 fflush(msgFile) ;
442 }
443 maxdomainsize=800;maxzeros=1000;maxsize=64;
444 /*maxdomainsize=neqns/100;*/
445 /*frontETree = orderViaMMD(graph, seed, msglvl, msgFile) ;*/
446 /*frontETree = orderViaND(graph,maxdomainsize,seed,msglvl,msgFile); */
447 /*frontETree = orderViaMS(graph,maxdomainsize,seed,msglvl,msgFile);*/
448 frontETree=orderViaBestOfNDandMS(graph,maxdomainsize,maxzeros,
449 maxsize,seed,msglvl,msgFile);
450 if ( msglvl > 1 ) {
451 fprintf(msgFile, "\n\n front tree from ordering") ;
452 ETree_writeForHumanEye(frontETree, msgFile) ;
453 fflush(msgFile) ;
454 }
455 /*--------------------------------------------------------------------*/
456 /*
457 -----------------------------------------------------
458 STEP 3: get the permutation, permute the matrix and
459 front tree and get the symbolic factorization
460 -----------------------------------------------------
461 */
462 oldToNewIV = ETree_oldToNewVtxPerm(frontETree) ;
463 oldToNew = IV_entries(oldToNewIV) ;
464 newToOldIV = ETree_newToOldVtxPerm(frontETree) ;
465 ETree_permuteVertices(frontETree, oldToNewIV) ;
466 InpMtx_permute(mtxA, oldToNew, oldToNew) ;
467 /* InpMtx_mapToUpperTriangle(mtxA) ;*/
468 InpMtx_changeCoordType(mtxA,INPMTX_BY_CHEVRONS);
469 InpMtx_changeStorageMode(mtxA,INPMTX_BY_VECTORS);
470 symbfacIVL = SymbFac_initFromInpMtx(frontETree, mtxA) ;
471 if ( msglvl > 1 ) {
472 fprintf(msgFile, "\n\n old-to-new permutation vector") ;
473 IV_writeForHumanEye(oldToNewIV, msgFile) ;
474 fprintf(msgFile, "\n\n new-to-old permutation vector") ;
475 IV_writeForHumanEye(newToOldIV, msgFile) ;
476 fprintf(msgFile, "\n\n front tree after permutation") ;
477 ETree_writeForHumanEye(frontETree, msgFile) ;
478 fprintf(msgFile, "\n\n input matrix after permutation") ;
479 InpMtx_writeForHumanEye(mtxA, msgFile) ;
480 fprintf(msgFile, "\n\n symbolic factorization") ;
481 IVL_writeForHumanEye(symbfacIVL, msgFile) ;
482 fflush(msgFile) ;
483 }
484 /*--------------------------------------------------------------------*/
485 /*
486 ------------------------------------------
487 STEP 4: initialize the front matrix object
488 ------------------------------------------
489 */
490 frontmtx = FrontMtx_new() ;
491 mtxmanager = SubMtxManager_new() ;
492 SubMtxManager_init(mtxmanager, NO_LOCK, 0) ;
493 FrontMtx_init(frontmtx, frontETree, symbfacIVL, type, symmetryflag,
494 FRONTMTX_DENSE_FRONTS, pivotingflag, NO_LOCK, 0, NULL,
495 mtxmanager, msglvl, msgFile) ;
496 /*--------------------------------------------------------------------*/
497 /*
498 -----------------------------------------
499 STEP 5: compute the numeric factorization
500 -----------------------------------------
501 */
502 chvmanager = ChvManager_new() ;
503 ChvManager_init(chvmanager, NO_LOCK, 1) ;
504 DVfill(10, cpus, 0.0) ;
505 IVfill(20, stats, 0) ;
506 rootchv = FrontMtx_factorInpMtx(frontmtx, mtxA, tau, 0.0, chvmanager,
507 &error,cpus, stats, msglvl, msgFile) ;
508 ChvManager_free(chvmanager) ;
509 if ( msglvl > 1 ) {
510 fprintf(msgFile, "\n\n factor matrix") ;
511 FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
512 fflush(msgFile) ;
513 }
514 if ( rootchv != NULL ) {
515 fprintf(msgFile, "\n\n matrix found to be singular\n") ;
516 exit(-1) ;
517 }
518 if(error>=0){
519 fprintf(msgFile,"\n\nerror encountered at front %d",error);
520 exit(-1);
521 }
522 /*--------------------------------------------------------------------*/
523 /*
524 --------------------------------------
525 STEP 6: post-process the factorization
526 --------------------------------------
527 */
528 FrontMtx_postProcess(frontmtx, msglvl, msgFile) ;
529 if ( msglvl > 1 ) {
530 fprintf(msgFile, "\n\n factor matrix after post-processing") ;
531 FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
532 fflush(msgFile) ;
533 }
534
535 /* reinitialize nodempc */
536
537 *mpcfree=1;
538 for(j=0;j<*nmpc;j++){
539 ipompc[j]=0;}
540
541 /* filling the RHS */
542
543 jrhs=0;
544 nrhs=1;
545 mtxB=DenseMtx_new();
546 mtxX=DenseMtx_new();
547
548 for(i=nindep;i>0;i--){
549 idof=indepdof[i-1];
550 if(ipointer[idof]>0){
551
552 /* new RHS column */
553
554 DenseMtx_init(mtxB, type, 0, 0, neqns, nrhs, 1, neqns) ;
555 DenseMtx_zero(mtxB) ;
556
557 index=ipointer[idof];
558 while(1){
559 irow=icoef[2*index-2]-1;
560 coef=xcoef[index-1];
561 DenseMtx_setRealEntry(mtxB,irow,jrhs,coef);
562 index=icoef[2*index-1];
563 if(index==0) break;
564 }
565
566 if ( msglvl > 1 ) {
567 fprintf(msgFile, "\n\n rhs matrix in original ordering") ;
568 DenseMtx_writeForHumanEye(mtxB, msgFile) ;
569 fflush(msgFile) ;
570 }
571
572 /*--------------------------------------------------------------------*/
573 /*
574 ---------------------------------------------------------
575 STEP 8: permute the right hand side into the new ordering
576 ---------------------------------------------------------
577 */
578 DenseMtx_permuteRows(mtxB, oldToNewIV) ;
579 if ( msglvl > 1 ) {
580 fprintf(msgFile, "\n\n right hand side matrix in new ordering") ;
581 DenseMtx_writeForHumanEye(mtxB, msgFile) ;
582 fflush(msgFile) ;
583 }
584 /*--------------------------------------------------------------------*/
585 /*
586 -------------------------------
587 STEP 9: solve the linear system
588 -------------------------------
589 */
590 DenseMtx_init(mtxX, type, 0, 0, neqns, nrhs, 1, neqns) ;
591 DenseMtx_zero(mtxX) ;
592 FrontMtx_solve(frontmtx, mtxX, mtxB, mtxmanager,cpus, msglvl, msgFile) ;
593 if ( msglvl > 1 ) {
594 fprintf(msgFile, "\n\n solution matrix in new ordering") ;
595 DenseMtx_writeForHumanEye(mtxX, msgFile) ;
596 fflush(msgFile) ;
597 }
598 /*--------------------------------------------------------------------*/
599 /*
600 --------------------------------------------------------
601 STEP 10: permute the solution into the original ordering
602 --------------------------------------------------------
603 */
604 DenseMtx_permuteRows(mtxX, newToOldIV) ;
605 if ( msglvl > 1 ) {
606 fprintf(msgFile, "\n\n solution matrix in original ordering") ;
607 DenseMtx_writeForHumanEye(mtxX, msgFile) ;
608 fflush(msgFile) ;
609 }
610
611
612 for(j=0;j<*nmpc;j++){
613 b=DenseMtx_entries(mtxX)[j];
614 if(fabs(b)>1.e-10){
615 nodempc[3**mpcfree-1]=ipompc[j];
616 node=(ITG)((idof+8)/8);
617 idir=idof+1-8*(node-1);
618 nodempc[3**mpcfree-3]=node;
619 nodempc[3**mpcfree-2]=idir;
620 coefmpc[*mpcfree-1]=b;
621 ipompc[j]=(*mpcfree)++;
622 if(*mpcfree>*memmpc_){
623 *memmpc_=(ITG)(1.1**memmpc_);
624 RENEW(nodempc,ITG,3**memmpc_);
625 RENEW(coefmpc,double,*memmpc_);
626 }
627 }
628 }
629 }
630 }
631 /*--------------------------------------------------------------------*/
632 /*
633 -----------
634 free memory
635 -----------
636 */
637 FrontMtx_free(frontmtx) ;
638 DenseMtx_free(mtxB) ;
639 DenseMtx_free(mtxX) ;
640 IV_free(newToOldIV) ;
641 IV_free(oldToNewIV) ;
642 InpMtx_free(mtxA) ;
643 ETree_free(frontETree) ;
644 IVL_free(symbfacIVL) ;
645 SubMtxManager_free(mtxmanager) ;
646 Graph_free(graph) ;
647
648 /* diagonal terms */
649
650 for(i=0;i<*nmpc;i++){
651 j=ilmpc[i]-1;
652 idof=ikmpc[i];
653 node=(ITG)((idof+7)/8);
654 idir=idof-8*(node-1);
655 nodempc[3**mpcfree-1]=ipompc[j];
656 nodempc[3**mpcfree-3]=node;
657 nodempc[3**mpcfree-2]=idir;
658 coefmpc[*mpcfree-1]=1.;
659 ipompc[j]=(*mpcfree)++;
660 if(*mpcfree>*memmpc_){
661 *memmpc_=(ITG)(1.1**memmpc_);
662 RENEW(nodempc,ITG,3**memmpc_);
663 RENEW(coefmpc,double,*memmpc_);
664 }
665 }
666
667 SFREE(ipointer);SFREE(indepdof);SFREE(icoef);SFREE(xcoef);
668
669 fclose(msgFile);
670
671 }
672 #endif
673
674 /* determining the effective size of nodempc and coefmpc for
675 the reallocation*/
676
677 *mpcend=0;
678 *mpcmult=0;
679 *maxlenmpc=0;
680 for(i=0;i<*nmpc;i++){
681 index=ipompc[i];
682 *mpcend=max(*mpcend,index);
683 nterm=1;
684 while(1){
685 index=nodempc[3*index-1];
686 if(index==0){
687 *mpcmult+=nterm*(nterm-1);
688 *maxlenmpc=max(*maxlenmpc,nterm);
689 break;
690 }
691 *mpcend=max(*mpcend,index);
692 nterm++;
693 }
694 }
695
696 SFREE(jmpc);
697
698 *nodempcp=nodempc;
699 *coefmpcp=coefmpc;
700
701 /* for(i=0;i<*nmpc;i++){
702 j=i+1;
703 FORTRAN(writempc,(ipompc,nodempc,coefmpc,labmpc,&j));
704 }*/
705
706 return;
707 }
708