root/src/cascade.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. cascade

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

/* [<][>][^][v][top][bottom][index][help] */