CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Macros | Functions
cascade.c File Reference
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "CalculiX.h"
Include dependency graph for cascade.c:

Go to the source code of this file.

Macros

#define min(a, b)   ((a) <= (b) ? (a) : (b))
 
#define max(a, b)   ((a) >= (b) ? (a) : (b))
 

Functions

void cascade (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, char *labmpc, ITG *nk, ITG *memmpc_, ITG *icascade, ITG *maxlenmpc, ITG *callfrommain, ITG *iperturb, ITG *ithermal)
 

Macro Definition Documentation

◆ max

#define max (   a,
  b 
)    ((a) >= (b) ? (a) : (b))

◆ min

#define min (   a,
  b 
)    ((a) <= (b) ? (a) : (b))

Function Documentation

◆ cascade()

void cascade ( 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,
char *  labmpc,
ITG *  nk,
ITG *  memmpc_,
ITG *  icascade,
ITG *  maxlenmpc,
ITG *  callfrommain,
ITG *  iperturb,
ITG *  ithermal 
)
38  {
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,ifluidmpc,
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  if(strcmp1(&labmpc[20*i],"FLUID")!=0){
113  printf("*ERROR in cascade: the DOF corresponding to \n node %" ITGFORMAT " in direction %" ITGFORMAT " is detected on the \n dependent side of a MPC and a SPC\n",
114  (ikmpc[i])/8+1,ikmpc[i]-8*((ikmpc[i])/8));
115  }else{
116  printf("*ERROR in cascade: the DOF corresponding to \n face %" ITGFORMAT " in direction %" ITGFORMAT " is detected on the \n dependent side of a MPC and a SPC\n",
117  (-ikmpc[i])/8+1,-ikmpc[i]-8*((-ikmpc[i])/8));
118  }
119  FORTRAN(stop,());
120  }
121  }
122  }
123 
124 /* check whether there are user mpc's: in user MPC's the
125  dependent DOF can change, however, the number of terms
126  cannot change */
127 
128  for(i=0;i<*nmpc;i++){
129 
130  /* linear mpc */
131 
132  /* because of the next line the size of field labmpc
133  has to be defined as 20*nmpc+1: without "+1" an
134  undefined field is accessed */
135 
136  if((strcmp1(&labmpc[20*i]," ")==0) ||
137  (strcmp1(&labmpc[20*i],"CYCLIC")==0) ||
138  (strcmp1(&labmpc[20*i],"SUBCYCLIC")==0)||
139  (strcmp1(&labmpc[20*i],"PRETENSION")==0)||
140  (strcmp1(&labmpc[20*i],"THERMALPRET")==0)||
141 // (strcmp1(&labmpc[20*i],"CONTACT")==0)||
142  (strcmp1(&labmpc[20*i],"FLUID")==0)||
143  (*iperturb<2)) jmpc[i]=0;
144 
145  /* nonlinear mpc */
146 
147  else if((strcmp1(&labmpc[20*i],"RIGID")==0) ||
148  (strcmp1(&labmpc[20*i],"KNOT")==0) ||
149  (strcmp1(&labmpc[20*i],"PLANE")==0) ||
150  (strcmp1(&labmpc[20*i],"BEAM")==0) ||
151  (strcmp1(&labmpc[20*i],"STRAIGHT")==0)) jmpc[i]=1;
152 
153  /* user mpc */
154 
155  else{
156  jmpc[i]=1;
157  if(*icascade==0) *icascade=1;
158  }
159  }
160 
161 /* decascading */
162 
163  ispooles=0;
164 
165  /* decascading using simple substitution */
166 
167  do{
168  ichange=0;
169  for(i=0;i<*nmpc;i++){
170 
171  if(strcmp1(&labmpc[20*i],"FLUID")!=0){
172  ifluidmpc=0;
173  }else{
174  ifluidmpc=1;
175  }
176 
177  if(jmpc[i]==1) nl=1;
178  else nl=0;
179  iexpand=0;
180  index=nodempc[3*ipompc[i]-1];
181  if(index==0) continue;
182  do{
183  if(ifluidmpc==0){
184  /* MPC on node */
185  idof=(nodempc[3*index-3]-1)*8+nodempc[3*index-2];
186  }else{
187  if(nodempc[3*index-3]>0){
188  /* MPC on face */
189  idof=-((nodempc[3*index-3]-1)*8+nodempc[3*index-2]);
190  }else{
191  /* MPC on node
192  SPC number: -nodempc[3*index-3]
193  node: nodeboun[-nodempc[3*index-3]-1] */
194 
195  idof=(nodeboun[-nodempc[3*index-3]-1]-1)*8+nodempc[3*index-2];
196  }
197  }
198 
199  FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
200  if((id>0)&&(ikmpc[id-1]==idof)){
201 
202  /* a term on the independent side of the MPC is
203  detected as dependent node in another MPC */
204 
205  indexold=nodempc[3*index-1];
206  coef=coefmpc[index-1];
207  mpc=ilmpc[id-1];
208 
209  /* no expansion if there is a dependence of a
210  nonlinear MPC on another linear or nonlinear MPC
211  and the call is from main */
212 
213  if((jmpc[mpc-1]==1)||(nl==1)){
214  *icascade=2;
215  if(idepend==0){
216  printf("*INFO in cascade: linear MPCs and\n");
217  printf(" nonlinear MPCs depend on each other\n");
218  printf(" common node: %" ITGFORMAT " in direction %" ITGFORMAT "\n\n",nodempc[3*index-3],nodempc[3*index-2]);
219  idepend=1;}
220  if(*callfrommain==1){
221  index=nodempc[3*index-1];
222  if(index!=0) continue;
223  else break;}
224  }
225 
226 /* printf("*INFO in cascade: DOF %" ITGFORMAT " of node %" ITGFORMAT " is expanded\n",
227  nodempc[3*index-2],nodempc[3*index-3]);*/
228 
229  /* collecting terms corresponding to the same DOF */
230 
231  index1=ipompc[i];
232  do{
233  index2old=index1;
234  index2=nodempc[3*index1-1];
235  if(index2==0) break;
236  do{
237  if((nodempc[3*index1-3]==nodempc[3*index2-3])&&
238  (nodempc[3*index1-2]==nodempc[3*index2-2])){
239  coefmpc[index1-1]+=coefmpc[index2-1];
240  nodempc[3*index2old-1]=nodempc[3*index2-1];
241  nodempc[3*index2-1]=*mpcfree;
242  *mpcfree=index2;
243  index2=nodempc[3*index2old-1];
244  if(index2==0) break;
245  }
246  else{
247  index2old=index2;
248  index2=nodempc[3*index2-1];
249  if(index2==0) break;
250  }
251  }while(1);
252  index1=nodempc[3*index1-1];
253  if(index1==0) break;
254  }while(1);
255 
256  /* check for zero coefficients on the dependent side */
257 
258  index1=ipompc[i];
259  if(fabs(coefmpc[index1-1])<1.e-10){
260  printf("*ERROR in cascade: zero coefficient on the\n");
261  printf(" dependent side of an equation\n");
262  printf(" dependent node: %" ITGFORMAT "",nodempc[3*index1-3]);
263  printf(" direction: %" ITGFORMAT "",nodempc[3*index1-2]);
264  FORTRAN(stop,());
265  }
266 
267  ichange=1;iexpand=1;
268  if((strcmp1(&labmpc[20*i]," ")==0)&&
269  (strcmp1(&labmpc[20*(mpc-1)],"CYCLIC")==0))
270  strcpy1(&labmpc[20*i],"SUBCYCLIC",9);
271  indexnew=ipompc[mpc-1];
272  coef=-coef/coefmpc[indexnew-1];
273  indexnew=nodempc[3*indexnew-1];
274  if(indexnew!=0){
275  do{
276  coefmpc[index-1]=coef*coefmpc[indexnew-1];
277  nodempc[3*index-3]=nodempc[3*indexnew-3];
278  nodempc[3*index-2]=nodempc[3*indexnew-2];
279  indexnew=nodempc[3*indexnew-1];
280  if(indexnew!=0){
281  nodempc[3*index-1]=*mpcfree;
282  index=*mpcfree;
283  *mpcfree=nodempc[3**mpcfree-1];
284  if(*mpcfree==0){
285  *mpcfree=*memmpc_+1;
286  nodempc[3*index-1]=*mpcfree;
287  *memmpc_=(ITG)(1.1**memmpc_);
288  printf("*INFO in cascade: reallocating nodempc; new size = %" ITGFORMAT "\n\n",*memmpc_);
289  RENEW(nodempc,ITG,3**memmpc_);
290  RENEW(coefmpc,double,*memmpc_);
291  for(j=*mpcfree;j<*memmpc_;j++){
292  nodempc[3*j-1]=j+1;
293  }
294  nodempc[3**memmpc_-1]=0;
295  }
296  continue;
297  }
298  else{
299  nodempc[3*index-1]=indexold;
300  break;
301  }
302  }while(1);
303  }else{
304  coefmpc[index-1]=0.;
305  }
306  break;
307  }
308  else{
309  index=nodempc[3*index-1];
310  if(index!=0) continue;
311  else break;
312  }
313  }while(1);
314  if(iexpand==0) continue;
315 
316  /* one term of the mpc was expanded
317  collecting terms corresponding to the same DOF */
318 
319  index1=ipompc[i];
320  do{
321  index2old=index1;
322  index2=nodempc[3*index1-1];
323  if(index2==0) break;
324  do{
325  if((nodempc[3*index1-3]==nodempc[3*index2-3])&&
326  (nodempc[3*index1-2]==nodempc[3*index2-2])){
327  coefmpc[index1-1]+=coefmpc[index2-1];
328  nodempc[3*index2old-1]=nodempc[3*index2-1];
329  nodempc[3*index2-1]=*mpcfree;
330  *mpcfree=index2;
331  index2=nodempc[3*index2old-1];
332  if(index2==0) break;
333  }
334  else{
335  index2old=index2;
336  index2=nodempc[3*index2-1];
337  if(index2==0) break;
338  }
339  }while(1);
340  index1=nodempc[3*index1-1];
341  if(index1==0) break;
342  }while(1);
343 
344  /* check for zero coefficients on the dependent and
345  independent side */
346 
347  index1=ipompc[i];
348  index1old=0;
349  do {
350  if(fabs(coefmpc[index1-1])<1.e-10){
351  if(index1old==0){
352  printf("*ERROR in cascade: zero coefficient on the\n");
353  printf(" dependent side of an equation\n");
354  printf(" dependent node: %" ITGFORMAT "",nodempc[3*index1-3]);
355  printf(" direction: %" ITGFORMAT "",nodempc[3*index1-2]);
356  FORTRAN(stop,());
357  }
358  else{
359  nodempc[3*index1old-1]=nodempc[3*index1-1];
360  nodempc[3*index1-1]=*mpcfree;
361  *mpcfree=index1;
362  index1=nodempc[3*index1old-1];
363  }
364  }
365  else{
366  index1old=index1;
367  index1=nodempc[3*index1-1];
368  }
369  if(index1==0) break;
370  }while(1);
371  }
372  if(ichange==0) break;
373  }while(1);
374 
375  /* decascading using spooles */
376 
377 #ifdef SPOOLES
378  if((*icascade==1)&&(ispooles==1)){
379  if ( (msgFile = fopen("spooles.out", "a")) == NULL ) {
380  fprintf(stderr, "\n fatal error in spooles.c"
381  "\n unable to open file spooles.out\n") ;
382  }
383  NNEW(ipointer,ITG,7**nk);
384  NNEW(indepdof,ITG,7**nk);
385  NNEW(icoef,ITG,2**memmpc_);
386  NNEW(xcoef,double,*memmpc_);
387  ifree=0;
388  nindep=0;
389 
390  for(i=*nmpc-1;i>-1;i--){
391  index=ipompc[i];
392  while(1){
393  idof=8*(nodempc[3*index-3]-1)+nodempc[3*index-2]-1;
394 
395 /* check whether idof is a independent dof which has not yet been
396  stored in indepdof */
397 
398  FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
399  if((id==0)||(ikmpc[id-1]!=idof)){
400  FORTRAN(nident,(indepdof,&idof,&nindep,&id));
401  if((id==0)||(indepdof[id-1]!=idof)){
402  for(j=nindep;j>id;j--){
403  indepdof[j]=indepdof[j-1];
404  }
405  indepdof[id]=idof;
406  nindep++;
407  }
408  }
409 
410  icoef[2*ifree]=i+1;
411  icoef[2*ifree+1]=ipointer[idof];
412  xcoef[ifree]=coefmpc[index-1];
413  ipointer[idof]=++ifree;
414  index=nodempc[3*index-1];
415  if(index==0) break;
416  }
417  }
418 
419 /* filling the left hand side */
420 
421  nent=*memmpc_;
422  neqns=*nmpc;
423  mtxA = InpMtx_new() ;
424  InpMtx_init(mtxA, INPMTX_BY_ROWS, type, nent, neqns) ;
425 
426  for(i=0;i<*nmpc;i++){
427  idof=ikmpc[i];
428  icolumn=ilmpc[i]-1;
429  if(strcmp1(&labmpc[20*icolumn],"RIGID")==0) icolnl=1;
430  else icolnl=0;
431  index=ipointer[idof-1];
432  while(1){
433  irow=icoef[2*index-2]-1;
434  if(irow!=icolumn){
435  if(strcmp1(&labmpc[20*irow],"RIGID")==0)irownl=1;
436  else irownl=0;
437  if((irownl==1)||(icolnl==1)){
438  *icascade=2;
439  InpMtx_free(mtxA);
440  printf("*ERROR in cascade: linear and nonlinear MPCs depend on each other");
441  FORTRAN(stop,());
442  }
443  }
444  if((strcmp1(&labmpc[20*irow]," ")==0)&&
445  (strcmp1(&labmpc[20*icolumn],"CYCLIC")==0)){
446  strcpy1(&labmpc[20*irow],"SUBCYCLIC",9);}
447  coef=xcoef[index-1];
448  InpMtx_inputRealEntry(mtxA,irow,icolumn,coef);
449  index=icoef[2*index-1];
450  if(index==0) break;
451  }
452  ipointer[idof-1]=0;
453  }
454 
455  InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ;
456  if ( msglvl > 1 ) {
457  fprintf(msgFile, "\n\n input matrix") ;
458  InpMtx_writeForHumanEye(mtxA, msgFile) ;
459  fflush(msgFile) ;
460  }
461 /*--------------------------------------------------------------------*/
462 /*
463  -------------------------------------------------
464  STEP 2 : find a low-fill ordering
465  (1) create the Graph object
466  (2) order the graph using multiple minimum degree
467  -------------------------------------------------
468 */
469  graph = Graph_new() ;
470  adjIVL = InpMtx_fullAdjacency(mtxA) ;
471  nedges = IVL_tsize(adjIVL) ;
472  Graph_init2(graph, 0, neqns, 0, nedges, neqns, nedges, adjIVL,
473  NULL, NULL) ;
474  if ( msglvl > 1 ) {
475  fprintf(msgFile, "\n\n graph of the input matrix") ;
476  Graph_writeForHumanEye(graph, msgFile) ;
477  fflush(msgFile) ;
478  }
479  maxdomainsize=800;maxzeros=1000;maxsize=64;
480  /*maxdomainsize=neqns/100;*/
481  /*frontETree = orderViaMMD(graph, seed, msglvl, msgFile) ;*/
482  /*frontETree = orderViaND(graph,maxdomainsize,seed,msglvl,msgFile); */
483  /*frontETree = orderViaMS(graph,maxdomainsize,seed,msglvl,msgFile);*/
484  frontETree=orderViaBestOfNDandMS(graph,maxdomainsize,maxzeros,
485  maxsize,seed,msglvl,msgFile);
486  if ( msglvl > 1 ) {
487  fprintf(msgFile, "\n\n front tree from ordering") ;
488  ETree_writeForHumanEye(frontETree, msgFile) ;
489  fflush(msgFile) ;
490  }
491 /*--------------------------------------------------------------------*/
492 /*
493  -----------------------------------------------------
494  STEP 3: get the permutation, permute the matrix and
495  front tree and get the symbolic factorization
496  -----------------------------------------------------
497 */
498  oldToNewIV = ETree_oldToNewVtxPerm(frontETree) ;
499  oldToNew = IV_entries(oldToNewIV) ;
500  newToOldIV = ETree_newToOldVtxPerm(frontETree) ;
501  ETree_permuteVertices(frontETree, oldToNewIV) ;
502  InpMtx_permute(mtxA, oldToNew, oldToNew) ;
503 /* InpMtx_mapToUpperTriangle(mtxA) ;*/
504  InpMtx_changeCoordType(mtxA,INPMTX_BY_CHEVRONS);
505  InpMtx_changeStorageMode(mtxA,INPMTX_BY_VECTORS);
506  symbfacIVL = SymbFac_initFromInpMtx(frontETree, mtxA) ;
507  if ( msglvl > 1 ) {
508  fprintf(msgFile, "\n\n old-to-new permutation vector") ;
509  IV_writeForHumanEye(oldToNewIV, msgFile) ;
510  fprintf(msgFile, "\n\n new-to-old permutation vector") ;
511  IV_writeForHumanEye(newToOldIV, msgFile) ;
512  fprintf(msgFile, "\n\n front tree after permutation") ;
513  ETree_writeForHumanEye(frontETree, msgFile) ;
514  fprintf(msgFile, "\n\n input matrix after permutation") ;
515  InpMtx_writeForHumanEye(mtxA, msgFile) ;
516  fprintf(msgFile, "\n\n symbolic factorization") ;
517  IVL_writeForHumanEye(symbfacIVL, msgFile) ;
518  fflush(msgFile) ;
519  }
520 /*--------------------------------------------------------------------*/
521 /*
522  ------------------------------------------
523  STEP 4: initialize the front matrix object
524  ------------------------------------------
525 */
526  frontmtx = FrontMtx_new() ;
527  mtxmanager = SubMtxManager_new() ;
528  SubMtxManager_init(mtxmanager, NO_LOCK, 0) ;
529  FrontMtx_init(frontmtx, frontETree, symbfacIVL, type, symmetryflag,
530  FRONTMTX_DENSE_FRONTS, pivotingflag, NO_LOCK, 0, NULL,
531  mtxmanager, msglvl, msgFile) ;
532 /*--------------------------------------------------------------------*/
533 /*
534  -----------------------------------------
535  STEP 5: compute the numeric factorization
536  -----------------------------------------
537 */
538  chvmanager = ChvManager_new() ;
539  ChvManager_init(chvmanager, NO_LOCK, 1) ;
540  DVfill(10, cpus, 0.0) ;
541  IVfill(20, stats, 0) ;
542  rootchv = FrontMtx_factorInpMtx(frontmtx, mtxA, tau, 0.0, chvmanager,
543  &error,cpus, stats, msglvl, msgFile) ;
544  ChvManager_free(chvmanager) ;
545  if ( msglvl > 1 ) {
546  fprintf(msgFile, "\n\n factor matrix") ;
547  FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
548  fflush(msgFile) ;
549  }
550  if ( rootchv != NULL ) {
551  fprintf(msgFile, "\n\n matrix found to be singular\n") ;
552  exit(-1) ;
553  }
554  if(error>=0){
555  fprintf(msgFile,"\n\nerror encountered at front %" ITGFORMAT "",error);
556  exit(-1);
557  }
558 /*--------------------------------------------------------------------*/
559 /*
560  --------------------------------------
561  STEP 6: post-process the factorization
562  --------------------------------------
563 */
564  FrontMtx_postProcess(frontmtx, msglvl, msgFile) ;
565  if ( msglvl > 1 ) {
566  fprintf(msgFile, "\n\n factor matrix after post-processing") ;
567  FrontMtx_writeForHumanEye(frontmtx, msgFile) ;
568  fflush(msgFile) ;
569  }
570 
571 /* reinitialize nodempc */
572 
573  *mpcfree=1;
574  for(j=0;j<*nmpc;j++){
575  ipompc[j]=0;}
576 
577 /* filling the RHS */
578 
579  jrhs=0;
580  nrhs=1;
581  mtxB=DenseMtx_new();
582  mtxX=DenseMtx_new();
583 
584  for(i=nindep;i>0;i--){
585  idof=indepdof[i-1];
586  if(ipointer[idof]>0){
587 
588 /* new RHS column */
589 
590  DenseMtx_init(mtxB, type, 0, 0, neqns, nrhs, 1, neqns) ;
591  DenseMtx_zero(mtxB) ;
592 
593  index=ipointer[idof];
594  while(1){
595  irow=icoef[2*index-2]-1;
596  coef=xcoef[index-1];
597  DenseMtx_setRealEntry(mtxB,irow,jrhs,coef);
598  index=icoef[2*index-1];
599  if(index==0) break;
600  }
601 
602  if ( msglvl > 1 ) {
603  fprintf(msgFile, "\n\n rhs matrix in original ordering") ;
604  DenseMtx_writeForHumanEye(mtxB, msgFile) ;
605  fflush(msgFile) ;
606  }
607 
608 /*--------------------------------------------------------------------*/
609 /*
610  ---------------------------------------------------------
611  STEP 8: permute the right hand side into the new ordering
612  ---------------------------------------------------------
613 */
614  DenseMtx_permuteRows(mtxB, oldToNewIV) ;
615  if ( msglvl > 1 ) {
616  fprintf(msgFile, "\n\n right hand side matrix in new ordering") ;
617  DenseMtx_writeForHumanEye(mtxB, msgFile) ;
618  fflush(msgFile) ;
619  }
620 /*--------------------------------------------------------------------*/
621 /*
622  -------------------------------
623  STEP 9: solve the linear system
624  -------------------------------
625 */
626  DenseMtx_init(mtxX, type, 0, 0, neqns, nrhs, 1, neqns) ;
627  DenseMtx_zero(mtxX) ;
628  FrontMtx_solve(frontmtx, mtxX, mtxB, mtxmanager,cpus, msglvl, msgFile) ;
629  if ( msglvl > 1 ) {
630  fprintf(msgFile, "\n\n solution matrix in new ordering") ;
631  DenseMtx_writeForHumanEye(mtxX, msgFile) ;
632  fflush(msgFile) ;
633  }
634 /*--------------------------------------------------------------------*/
635 /*
636  --------------------------------------------------------
637  STEP 10: permute the solution into the original ordering
638  --------------------------------------------------------
639 */
640  DenseMtx_permuteRows(mtxX, newToOldIV) ;
641  if ( msglvl > 1 ) {
642  fprintf(msgFile, "\n\n solution matrix in original ordering") ;
643  DenseMtx_writeForHumanEye(mtxX, msgFile) ;
644  fflush(msgFile) ;
645  }
646 
647 
648  for(j=0;j<*nmpc;j++){
649  b=DenseMtx_entries(mtxX)[j];
650  if(fabs(b)>1.e-10){
651  nodempc[3**mpcfree-1]=ipompc[j];
652  node=(ITG)((idof+8)/8);
653  idir=idof+1-8*(node-1);
654  nodempc[3**mpcfree-3]=node;
655  nodempc[3**mpcfree-2]=idir;
656  coefmpc[*mpcfree-1]=b;
657  ipompc[j]=(*mpcfree)++;
658  if(*mpcfree>*memmpc_){
659  *memmpc_=(ITG)(1.1**memmpc_);
660  RENEW(nodempc,ITG,3**memmpc_);
661  RENEW(coefmpc,double,*memmpc_);
662  }
663  }
664  }
665  }
666  }
667 /*--------------------------------------------------------------------*/
668 /*
669  -----------
670  free memory
671  -----------
672 */
673  FrontMtx_free(frontmtx) ;
674  DenseMtx_free(mtxB) ;
675  DenseMtx_free(mtxX) ;
676  IV_free(newToOldIV) ;
677  IV_free(oldToNewIV) ;
678  InpMtx_free(mtxA) ;
679  ETree_free(frontETree) ;
680  IVL_free(symbfacIVL) ;
681  SubMtxManager_free(mtxmanager) ;
682  Graph_free(graph) ;
683 
684 /* diagonal terms */
685 
686  for(i=0;i<*nmpc;i++){
687  j=ilmpc[i]-1;
688  idof=ikmpc[i];
689  node=(ITG)((idof+7)/8);
690  idir=idof-8*(node-1);
691  nodempc[3**mpcfree-1]=ipompc[j];
692  nodempc[3**mpcfree-3]=node;
693  nodempc[3**mpcfree-2]=idir;
694  coefmpc[*mpcfree-1]=1.;
695  ipompc[j]=(*mpcfree)++;
696  if(*mpcfree>*memmpc_){
697  *memmpc_=(ITG)(1.1**memmpc_);
698  RENEW(nodempc,ITG,3**memmpc_);
699  RENEW(coefmpc,double,*memmpc_);
700  }
701  }
702 
703  SFREE(ipointer);SFREE(indepdof);SFREE(icoef);SFREE(xcoef);
704 
705  fclose(msgFile);
706 
707  }
708 #endif
709 
710 /* determining the effective size of nodempc and coefmpc for
711  the reallocation*/
712 
713  *mpcend=0;
714  *maxlenmpc=0;
715  for(i=0;i<*nmpc;i++){
716  index=ipompc[i];
717  *mpcend=max(*mpcend,index);
718  nterm=1;
719  while(1){
720  index=nodempc[3*index-1];
721  if(index==0){
722  *maxlenmpc=max(*maxlenmpc,nterm);
723  break;
724  }
725  *mpcend=max(*mpcend,index);
726  nterm++;
727  }
728  }
729 
730  SFREE(jmpc);
731 
732  *nodempcp=nodempc;
733  *coefmpcp=coefmpc;
734 
735  /* for(i=0;i<*nmpc;i++){
736  j=i+1;
737  FORTRAN(writempc,(ipompc,nodempc,coefmpc,labmpc,&j));
738  }*/
739 
740  return;
741 }
ITGFORMAT
#define ITGFORMAT
Definition: CalculiX.h:52
max
#define max(a, b)
Definition: cascade.c:32
strcpy1
ITG strcpy1(char *s1, const char *s2, ITG length)
Definition: strcpy1.c:24
FORTRAN
void FORTRAN(actideacti,(char *set, ITG *nset, ITG *istartset, ITG *iendset, ITG *ialset, char *objectset, ITG *ipkon, ITG *ibject, ITG *ne))
strcmp1
ITG strcmp1(const char *s1, const char *s2)
Definition: strcmp1.c:24
stop
subroutine stop()
Definition: stop.f:20
tau
void tau(double *ad, double **aup, double *adb, double *aubp, double *sigma, double *b, ITG *icol, ITG **irowp, ITG *neq, ITG *nzs)
RENEW
#define RENEW(a, b, c)
Definition: CalculiX.h:40
SFREE
#define SFREE(a)
Definition: CalculiX.h:41
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
ITG
#define ITG
Definition: CalculiX.h:51
NNEW
#define NNEW(a, b, c)
Definition: CalculiX.h:39
Hosted by OpenAircraft.com, (Michigan UAV, LLC)