root/src/arpack.c

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

DEFINITIONS

This source file includes following definitions.
  1. arpack

   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 #ifdef ARPACK
  19 
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include "CalculiX.h"
  25 #ifdef SPOOLES
  26    #include "spooles.h"
  27 #endif
  28 #ifdef SGI
  29    #include "sgi.h"
  30 #endif
  31 #ifdef TAUCS
  32    #include "tau.h"
  33 #endif
  34 #ifdef MATRIXSTORAGE
  35    #include "matrixstorage.h"
  36 #endif
  37 #ifdef PARDISO
  38    #include "pardiso.h"
  39 #endif
  40 
  41 void arpack(double *co, ITG *nk, ITG **konp, ITG **ipkonp, char **lakonp,
  42              ITG *ne, 
  43              ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, 
  44              ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc, 
  45              ITG *nmpc, 
  46              ITG *nodeforc, ITG *ndirforc,double *xforc, ITG *nforc, 
  47              ITG *nelemload, char *sideload, double *xload,
  48              ITG *nload, ITG *nactdof, 
  49              ITG *icol, ITG *jq, ITG **irowp, ITG *neq, ITG *nzl, 
  50              ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, 
  51              ITG *ilboun,
  52              double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon,
  53              double *shcon, ITG *nshcon, double *cocon, ITG *ncocon,
  54              double *alcon, ITG *nalcon, double *alzero, ITG **ielmatp,
  55              ITG **ielorienp, ITG *norien, double *orab, ITG *ntmat_,
  56              double *t0, double *t1, double *t1old, 
  57              ITG *ithermal,double *prestr, ITG *iprestr, 
  58              double *vold,ITG *iperturb, double *sti, ITG *nzs,  
  59              ITG *kode, ITG *mei, double *fei,
  60              char *filab, double *eme,
  61              ITG *iexpl, double *plicon, ITG *nplicon, double *plkcon,
  62              ITG *nplkcon,
  63              double **xstatep, ITG *npmat_, char *matname, ITG *mi,
  64              ITG *ncmat_, ITG *nstate_, double **enerp, char *jobnamec,
  65              char *output, char *set, ITG *nset, ITG *istartset,
  66              ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
  67              char *prset, ITG *nener, ITG *isolver, double *trab, 
  68              ITG *inotr, ITG *ntrans, double *ttime, double *fmpc,
  69              char *cbody, ITG *ibody,double *xbody, ITG *nbody,
  70              double *thicke, ITG *nslavs, double *tietol, ITG *nkon,
  71              ITG *mpcinfo,ITG *ntie,ITG *istep,ITG *mcs,ITG *ics,
  72              char *tieset,double *cs,ITG *nintpoint,ITG *mortar,
  73              ITG *ifacecount,ITG **islavsurfp,double **pslavsurfp,
  74              double **clearinip,ITG *nmat,char *typeboun,
  75              ITG *ielprop,double *prop){
  76 
  77   /* calls the Arnoldi Package (ARPACK) */
  78   
  79   char bmat[2]="G", which[3]="LM", howmny[2]="A", fneig[132]="",
  80       description[13]="            ",*lakon=NULL;
  81 
  82   ITG *inum=NULL,k,ido,ldz,iparam[11],ipntr[14],lworkl,ngraph=1,im,
  83     info,rvec=1,*select=NULL,lfin,j,lint,iout,ielas=0,icmd=0,mt=mi[1]+1,
  84     iinc=1,nev,ncv,mxiter,jrow,*ipobody=NULL,inewton=0,ifreebody,
  85     mass[2]={1,1}, stiffness=1, buckling=0, rhsi=0, intscheme=0,noddiam=-1,
  86     coriolis=0,symmetryflag=0,inputformat=0,*ipneigh=NULL,*neigh=NULL,ne0,
  87     *integerglob=NULL,nasym=0,zero=0,irenewxstate,ncont=0,*itietri=NULL,
  88     *koncont=NULL,ismallsliding=0,*itiefac=NULL,*islavsurf=NULL,
  89     *islavnode=NULL,*imastnode=NULL,*nslavnode=NULL,*nmastnode=NULL,
  90     *imastop=NULL,*iponoels=NULL,*inoels=NULL,*ipe=NULL,*ime=NULL,
  91     mpcfree,memmpc_,icascade,maxlenmpc,nkon0,iit=-1,*irow=NULL,nherm=1,
  92     icfd=0,*inomat=NULL,*ipkon=NULL,*kon=NULL,*ielmat=NULL,*ielorien=NULL,
  93     *islavact=NULL,maxprevcontel,iex,nslavs_prev_step,
  94     iflagact=0,*islavsurfold=NULL;
  95 
  96   double *stn=NULL,*v=NULL,*resid=NULL,*z=NULL,*workd=NULL,
  97     *workl=NULL,*d=NULL,sigma=1,*temp_array=NULL,*pslavsurfold=NULL,
  98     *een=NULL,sum,cam[5],*f=NULL,*fn=NULL,qa[3],*fext=NULL,
  99     *epn=NULL,*fnr=NULL,*fni=NULL,*emn=NULL,*emeini=NULL,
 100     *xstateini=NULL,*xstiff=NULL,*stiini=NULL,*vini=NULL,freq,*stx=NULL,
 101     *enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,
 102     *physcon=NULL,*qfx=NULL,*qfn=NULL,tol,fmin,fmax,pi,*cgr=NULL,
 103     *xloadold=NULL,reltime=1.,*vr=NULL,*vi=NULL,*stnr=NULL,*stni=NULL,
 104     *vmax=NULL,*stnmax=NULL,*springarea=NULL,*eenmax=NULL,
 105     *doubleglob=NULL,*cg=NULL,*straight=NULL,*clearini=NULL,
 106     *xmastnor=NULL,*areaslav=NULL,*xnoels=NULL,theta=0.,
 107     *di=NULL,sigmai=0,*workev=NULL,*ener=NULL,*xstate=NULL,*dc=NULL,
 108     *au=NULL,*ad=NULL,*b=NULL,*aub=NULL,*adb=NULL,*pslavsurf=NULL,
 109     *pmastsurf=NULL,*cdn=NULL,
 110     *cdnr=NULL,*cdni=NULL;
 111 
 112   FILE *f1;
 113 
 114   /* dummy arguments for the results call */
 115 
 116   double *veold=NULL,*accold=NULL,bet,gam,dtime,time=1.;
 117 
 118 #ifdef SGI
 119   ITG token;
 120 #endif
 121 
 122   irow=*irowp;ener=*enerp;xstate=*xstatep;ipkon=*ipkonp;lakon=*lakonp;
 123   kon=*konp;ielmat=*ielmatp;ielorien=*ielorienp;
 124 
 125   islavsurf=*islavsurfp;pslavsurf=*pslavsurfp;clearini=*clearinip;
 126   
 127   if(*mortar==0){
 128       maxprevcontel=*nslavs;
 129   }else if(*mortar==1){
 130       maxprevcontel=*nintpoint;
 131   }
 132   nslavs_prev_step=*nslavs;
 133 
 134   if((strcmp1(&filab[870],"PU  ")==0)||
 135      (strcmp1(&filab[1479],"PHS ")==0)||
 136      (strcmp1(&filab[1566],"MAXU")==0)||
 137      (strcmp1(&filab[1653],"MAXS")==0)){
 138       printf("*WARNING in arpack: PU, PHS, MAXU or MAX was selected in a frequency calculation without cyclic symmetry;\n this is not correct; output request is removed;\n");
 139       strcpy1(&filab[870],"    ",4);
 140       strcpy1(&filab[1479],"    ",4);
 141       strcpy1(&filab[1566],"    ",4);
 142       strcpy1(&filab[1653],"    ",4);
 143   }
 144 
 145   /* copying the frequency parameters */
 146 
 147   pi=4.*atan(1.);
 148 
 149   nev=mei[0];
 150   ncv=mei[1];
 151   mxiter=mei[2];
 152   tol=fei[0];
 153   fmin=2*pi*fei[1];
 154   fmax=2*pi*fei[2];
 155   
 156   /* assigning the body forces to the elements */ 
 157 
 158   if(*nbody>0){
 159       ifreebody=*ne+1;
 160       NNEW(ipobody,ITG,2*ifreebody**nbody);
 161       for(k=1;k<=*nbody;k++){
 162           FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
 163                              iendset,ialset,&inewton,nset,&ifreebody,&k));
 164           RENEW(ipobody,ITG,2*(*ne+ifreebody));
 165       }
 166       RENEW(ipobody,ITG,2*(ifreebody-1));
 167       if(inewton==1){
 168           printf("*ERROR in arpackcs: generalized gravity loading is not allowed in frequency calculations");
 169       FORTRAN(stop,());
 170       }
 171   }
 172 
 173   ne0=*ne;nkon0=*nkon;
 174 
 175   /* contact conditions */
 176   
 177   if(*iperturb!=0){
 178 
 179       memmpc_=mpcinfo[0];mpcfree=mpcinfo[1];icascade=mpcinfo[2];
 180       maxlenmpc=mpcinfo[3];
 181 
 182       if(*nslavs==0){irenewxstate=1;}else{irenewxstate=0;}
 183       inicont(nk,&ncont,ntie,tieset,nset,set,istartset,iendset,ialset,&itietri,
 184           lakon,ipkon,kon,&koncont,nslavs,tietol,&ismallsliding,&itiefac,
 185           &islavsurf,&islavnode,&imastnode,&nslavnode,&nmastnode,
 186           mortar,&imastop,nkon,&iponoels,&inoels,&ipe,&ime,ne,ifacecount,
 187           nmpc,&mpcfree,&memmpc_,
 188           &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
 189           iperturb,ikboun,nboun,co,istep,&xnoels);
 190 
 191       if(ncont!=0){
 192 
 193           NNEW(cg,double,3*ncont);
 194           NNEW(straight,double,16*ncont);
 195           
 196           /* 11 instead of 10: last position is reserved for the
 197              local contact spring element number; needed as
 198              pointer into springarea */
 199           
 200           if(*mortar==0){
 201               RENEW(kon,ITG,*nkon+11**nslavs);
 202               NNEW(springarea,double,2**nslavs);
 203               if(*nener==1){
 204                   RENEW(ener,double,mi[0]*(*ne+*nslavs)*2);
 205               }
 206               RENEW(ipkon,ITG,*ne+*nslavs);
 207               RENEW(lakon,char,8*(*ne+*nslavs));
 208               
 209               if(*norien>0){
 210                   RENEW(ielorien,ITG,mi[2]*(*ne+*nslavs));
 211                   for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielorien[k]=0;
 212               }
 213               
 214               RENEW(ielmat,ITG,mi[2]*(*ne+*nslavs));
 215               for(k=mi[2]**ne;k<mi[2]*(*ne+*nslavs);k++) ielmat[k]=1;
 216               
 217               if((maxprevcontel==0)&&(*nslavs!=0)){
 218                   RENEW(xstate,double,*nstate_*mi[0]*(*ne+*nslavs));
 219                   for(k=*nstate_*mi[0]**ne;k<*nstate_*mi[0]*(*ne+*nslavs);k++){
 220                       xstate[k]=0.;
 221                   }
 222               }
 223               maxprevcontel=*nslavs;
 224               
 225               NNEW(areaslav,double,*ifacecount);
 226               NNEW(xmastnor,double,3*nmastnode[*ntie]);
 227           }else if(*mortar==1){
 228               NNEW(islavact,ITG,nslavnode[*ntie]);
 229               DMEMSET(islavact,0,nslavnode[*ntie],1);
 230               if((*istep==1)||(nslavs_prev_step==0)) NNEW(clearini,double,3*9**ifacecount);
 231               NNEW(xmastnor,double,3*nmastnode[*ntie]);
 232 
 233               if(*nstate_!=0){
 234                   if(maxprevcontel!=0){
 235                       NNEW(islavsurfold,ITG,2**ifacecount+2);
 236                       NNEW(pslavsurfold,double,3**nintpoint);
 237                       memcpy(&islavsurfold[0],&islavsurf[0],
 238                              sizeof(ITG)*(2**ifacecount+2));
 239                       memcpy(&pslavsurfold[0],&pslavsurf[0],
 240                              sizeof(double)*(3**nintpoint));
 241                   }
 242               }
 243 
 244               *nintpoint=0;
 245               
 246               precontact(&ncont,ntie,tieset,nset,set,istartset,
 247                          iendset,ialset,itietri,lakon,ipkon,kon,koncont,ne,
 248                          cg,straight,co,vold,istep,&iinc,&iit,itiefac,
 249                          islavsurf,islavnode,imastnode,nslavnode,nmastnode,
 250                          imastop,mi,ipe,ime,tietol,&iflagact,
 251                          nintpoint,&pslavsurf,xmastnor,cs,mcs,ics,clearini,
 252                          nslavs);
 253               
 254               /* changing the dimension of element-related fields */
 255               
 256               RENEW(kon,ITG,*nkon+22**nintpoint);
 257               RENEW(springarea,double,2**nintpoint);
 258               RENEW(pmastsurf,double,6**nintpoint);
 259               
 260               if(*nener==1){
 261                   RENEW(ener,double,mi[0]*(*ne+*nintpoint)*2);
 262               }
 263               RENEW(ipkon,ITG,*ne+*nintpoint);
 264               RENEW(lakon,char,8*(*ne+*nintpoint));
 265               
 266               if(*norien>0){
 267                   RENEW(ielorien,ITG,mi[2]*(*ne+*nintpoint));
 268                   for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielorien[k]=0;
 269               }
 270               RENEW(ielmat,ITG,mi[2]*(*ne+*nintpoint));
 271               for(k=mi[2]**ne;k<mi[2]*(*ne+*nintpoint);k++) ielmat[k]=1;
 272           }
 273           
 274           /* generating contact spring elements */
 275 
 276           contact(&ncont,ntie,tieset,nset,set,istartset,iendset,
 277              ialset,itietri,lakon,ipkon,kon,koncont,ne,cg,straight,nkon,
 278              co,vold,ielmat,cs,elcon,istep,&iinc,&iit,ncmat_,ntmat_,
 279              &ne0,vini,nmethod,nmpc,&mpcfree,&memmpc_,
 280              &ipompc,&labmpc,&ikmpc,&ilmpc,&fmpc,&nodempc,&coefmpc,
 281              iperturb,ikboun,nboun,mi,imastop,nslavnode,islavnode,islavsurf,
 282              itiefac,areaslav,iponoels,inoels,springarea,tietol,&reltime,
 283              imastnode,nmastnode,xmastnor,filab,mcs,ics,&nasym,
 284              xnoels,mortar,pslavsurf,pmastsurf,clearini,&theta);
 285           
 286           printf("number of contact spring elements=%" ITGFORMAT "\n\n",*ne-ne0);
 287                   
 288           /* interpolating the state variables */
 289 
 290           if(*mortar==1){
 291               if(*nstate_!=0){
 292                   if(maxprevcontel!=0){
 293                       RENEW(xstateini,double,
 294                             *nstate_*mi[0]*(ne0+maxprevcontel));
 295                       for(k=*nstate_*mi[0]*ne0;
 296                           k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
 297                           xstateini[k]=xstate[k];
 298                       }
 299                   }
 300                   
 301                   RENEW(xstate,double,*nstate_*mi[0]*(ne0+*nintpoint));
 302                   for(k=*nstate_*mi[0]*ne0;k<*nstate_*mi[0]*(ne0+*nintpoint);k++){
 303                       xstate[k]=0.;
 304                   }
 305                   
 306                   if((*nintpoint>0)&&(maxprevcontel>0)){
 307                       iex=2;
 308                       
 309                       /* interpolation of xstate */
 310                       
 311                       FORTRAN(interpolatestate,(ne,ipkon,kon,lakon,
 312                                &ne0,mi,xstate,pslavsurf,nstate_,
 313                                xstateini,islavsurf,islavsurfold,
 314                                pslavsurfold));
 315                       
 316                   }
 317 
 318                   if(maxprevcontel!=0){
 319                       SFREE(islavsurfold);SFREE(pslavsurfold);
 320                   }
 321 
 322                   maxprevcontel=*nintpoint;
 323                   
 324                   RENEW(xstateini,double,*nstate_*mi[0]*(ne0+*nintpoint));
 325                   for(k=0;k<*nstate_*mi[0]*(ne0+*nintpoint);++k){
 326                       xstateini[k]=xstate[k];
 327                   }
 328               }
 329           }
 330           
 331           /* determining the structure of the stiffness/mass matrix */
 332           
 333           remastructar(ipompc,&coefmpc,&nodempc,nmpc,
 334                  &mpcfree,nodeboun,ndirboun,nboun,ikmpc,ilmpc,ikboun,ilboun,
 335                  labmpc,nk,&memmpc_,&icascade,&maxlenmpc,
 336                  kon,ipkon,lakon,ne,nactdof,icol,jq,&irow,isolver,
 337                  neq,nzs,nmethod,ithermal,iperturb,mass,mi,ics,cs,
 338                  mcs,mortar,typeboun);
 339       }
 340   }
 341 
 342   /* field for initial values of state variables (needed if
 343      previous static step was viscoplastic and for contact */
 344 
 345   if((*nstate_!=0)&&((*mortar==0)||(ncont==0))){
 346       NNEW(xstateini,double,*nstate_*mi[0]*(ne0+*nslavs));
 347       for(k=0;k<*nstate_*mi[0]*(ne0+*nslavs);++k){
 348       xstateini[k]=xstate[k];
 349     }
 350   }
 351 
 352   /* determining the internal forces and the stiffness coefficients */
 353 
 354   NNEW(f,double,neq[1]);
 355 
 356   /* allocating a field for the stiffness matrix */
 357 
 358   NNEW(xstiff,double,(long long)27*mi[0]**ne);
 359 
 360   iout=-1;
 361   NNEW(v,double,mt**nk);
 362   memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
 363   NNEW(fn,double,mt**nk);
 364   NNEW(stx,double,6*mi[0]**ne);
 365   if(*ithermal>1){
 366       NNEW(qfx,double,3*mi[0]**ne);
 367   }
 368   NNEW(inum,ITG,*nk);
 369   if(*iperturb==0){
 370      results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 371              elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 372              ielorien,norien,orab,ntmat_,t0,t0,ithermal,
 373              prestr,iprestr,filab,eme,emn,een,iperturb,
 374              f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 375              ndirboun,xboun,nboun,ipompc,
 376              nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
 377              &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 378              xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 379              &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
 380              emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 381              iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 382              fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
 383              &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
 384              sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 385              mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 386              islavsurf,ielprop,prop);
 387   }else{
 388      results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
 389              elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 390              ielorien,norien,orab,ntmat_,t0,t1old,ithermal,
 391              prestr,iprestr,filab,eme,emn,een,iperturb,
 392              f,fn,nactdof,&iout,qa,vold,b,nodeboun,
 393              ndirboun,xboun,nboun,ipompc,
 394              nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
 395              &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 396              xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 397              &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
 398              emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 399              iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 400              fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
 401              &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
 402              sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 403              mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 404              islavsurf,ielprop,prop);
 405   }
 406   SFREE(f);SFREE(v);SFREE(fn);SFREE(stx);if(*ithermal>1)SFREE(qfx);SFREE(inum);
 407   iout=1;
 408 
 409   /* for the frequency analysis linear strain and elastic properties
 410      are used */
 411 
 412   iperturb[1]=0;ielas=1;
 413 
 414   /* filling in the matrix */
 415 
 416   NNEW(ad,double,neq[1]);
 417   NNEW(au,double,nzs[2]);
 418 
 419   NNEW(adb,double,neq[1]);
 420   NNEW(aub,double,nzs[1]);
 421 
 422   NNEW(fext,double,neq[1]);
 423 
 424   if(*iperturb==0){
 425     FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 426               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 427               nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 428               ad,au,fext,nactdof,icol,jq,irow,neq,nzl,nmethod,
 429               ikmpc,ilmpc,ikboun,ilboun,
 430               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 431               ielorien,norien,orab,ntmat_,
 432               t0,t0,ithermal,prestr,iprestr,vold,iperturb,sti,
 433               nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 434               xstiff,npmat_,&dtime,matname,mi,
 435               ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
 436               physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
 437               &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
 438               xstateini,xstate,thicke,integerglob,doubleglob,
 439               tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
 440               pmastsurf,mortar,clearini,ielprop,prop));
 441   }
 442   else{
 443       FORTRAN(mafillsm,(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xboun,nboun,
 444               ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 445               nforc,nelemload,sideload,xload,nload,xbody,ipobody,nbody,cgr,
 446               ad,au,fext,nactdof,icol,jq,irow,neq,nzl,nmethod,
 447               ikmpc,ilmpc,ikboun,ilboun,
 448               elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
 449               ielorien,norien,orab,ntmat_,
 450               t0,t1old,ithermal,prestr,iprestr,vold,iperturb,sti,
 451               nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 452               xstiff,npmat_,&dtime,matname,mi,
 453               ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
 454               physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
 455               &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
 456               xstateini,xstate,thicke,integerglob,doubleglob,
 457               tieset,istartset,iendset,ialset,ntie,&nasym,pslavsurf,
 458               pmastsurf,mortar,clearini,ielprop,prop));
 459 
 460       if(nasym==1){
 461           RENEW(au,double,nzs[2]+nzs[1]);
 462           RENEW(aub,double,nzs[2]+nzs[1]);
 463           symmetryflag=2;
 464           inputformat=1;
 465           
 466           FORTRAN(mafillsmas,(co,nk,kon,ipkon,lakon,ne,nodeboun,
 467                   ndirboun,xboun,nboun,
 468                   ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
 469                   nforc,nelemload,sideload,xload,nload,xbody,ipobody,
 470                   nbody,cgr,ad,au,fext,nactdof,icol,jq,irow,neq,nzl,
 471                   nmethod,ikmpc,ilmpc,ikboun,ilboun,
 472                   elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
 473                   ielmat,ielorien,norien,orab,ntmat_,
 474                   t0,t1old,ithermal,prestr,iprestr,vold,iperturb,sti,
 475                   nzs,stx,adb,aub,iexpl,plicon,nplicon,plkcon,nplkcon,
 476                   xstiff,npmat_,&dtime,matname,mi,
 477                   ncmat_,mass,&stiffness,&buckling,&rhsi,&intscheme,
 478                   physcon,shcon,nshcon,cocon,ncocon,ttime,&time,istep,&iinc,
 479                   &coriolis,ibody,xloadold,&reltime,veold,springarea,nstate_,
 480                   xstateini,xstate,thicke,
 481                   integerglob,doubleglob,tieset,istartset,iendset,
 482                   ialset,ntie,&nasym,pslavsurf,pmastsurf,mortar,clearini,
 483                   ielprop,prop));
 484       }
 485   }
 486 
 487   SFREE(fext);
 488 
 489   if(*nmethod==0){
 490 
 491     /* error occurred in mafill: storing the geometry in frd format */
 492 
 493     ++*kode;
 494     NNEW(inum,ITG,*nk);for(k=0;k<*nk;k++) inum[k]=1;
 495     if(strcmp1(&filab[1044],"ZZS")==0){
 496         NNEW(neigh,ITG,40**ne);
 497         NNEW(ipneigh,ITG,*nk);
 498     }
 499 
 500     frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
 501             kode,filab,een,t1,fn,&time,epn,ielmat,matname,enern,xstaten,
 502             nstate_,istep,&iinc,ithermal,qfn,&j,&noddiam,trab,inotr,
 503             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 504             mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 505             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 506             thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
 507     
 508     if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
 509     SFREE(inum);FORTRAN(stop,());
 510 
 511   }
 512 
 513   /* LU decomposition of the left hand matrix */
 514 
 515   if(nasym==1){sigma=0.;}else{sigma=1.;}
 516 
 517   if(*isolver==0){
 518 #ifdef SPOOLES
 519     spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
 520                    &symmetryflag,&inputformat,&nzs[2]);
 521 #else
 522     printf("*ERROR in arpack: the SPOOLES library is not linked\n\n");
 523     FORTRAN(stop,());
 524 #endif
 525   }
 526   else if(*isolver==4){
 527 #ifdef SGI
 528     token=1;
 529     sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
 530 #else
 531     printf("*ERROR in arpack: the SGI library is not linked\n\n");
 532     FORTRAN(stop,());
 533 #endif
 534   }
 535   else if(*isolver==5){
 536 #ifdef TAUCS
 537     tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
 538 #else
 539     printf("*ERROR in arpack: the TAUCS library is not linked\n\n");
 540     FORTRAN(stop,());
 541 #endif
 542   }
 543   else if(*isolver==6){
 544 #ifdef MATRIXSTORAGE
 545     matrixstorage(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1],
 546                   ntrans,inotr,trab,co,nk,nactdof,jobnamec,mi,ipkon,
 547                   lakon,kon,ne,mei,nboun,nmpc,cs,mcs);
 548 #else
 549     printf("*ERROR in arpack: the MATRIXSTORAGE library is not linked\n\n");
 550     FORTRAN(stop,());
 551 #endif
 552   }
 553   else if(*isolver==7){
 554 #ifdef PARDISO
 555     pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
 556                    &symmetryflag,&inputformat,jq,&nzs[2]);
 557 #else
 558     printf("*ERROR in arpack: the PARDISO library is not linked\n\n");
 559     FORTRAN(stop,());
 560 #endif
 561   }
 562 
 563 /* calculating the eigenvalues and eigenmodes */
 564   
 565   printf(" Calculating the eigenvalues and the eigenmodes\n\n");
 566 
 567   ido=0;
 568   ldz=neq[1];
 569   iparam[0]=1;
 570   iparam[2]=mxiter;
 571   iparam[3]=1;
 572   iparam[6]=3;
 573 
 574   info=0;
 575 
 576   NNEW(resid,double,neq[1]);
 577   NNEW(z,double,(long long)ncv*neq[1]);
 578   NNEW(workd,double,3*neq[1]);
 579 
 580   if(nasym==1){
 581       lworkl=3*ncv*(2+ncv);
 582       NNEW(workl,double,lworkl);
 583       FORTRAN(dnaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,
 584           workl,&lworkl,&info));
 585   }else{
 586       lworkl=ncv*(8+ncv);
 587       NNEW(workl,double,lworkl);
 588       FORTRAN(dsaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,
 589           workl,&lworkl,&info));
 590   }
 591 
 592   NNEW(temp_array,double,neq[1]);
 593 
 594   while((ido==-1)||(ido==1)||(ido==2)){
 595     if(ido==-1){
 596         if(nasym==1){
 597             FORTRAN(opas,(&neq[1],&workd[ipntr[0]-1],temp_array,adb,aub,jq,irow,nzs));
 598         }else{
 599             FORTRAN(op,(&neq[1],&workd[ipntr[0]-1],temp_array,adb,aub,jq,irow));
 600         }
 601     }
 602     if((ido==-1)||(ido==1)){
 603 
 604       /* solve the linear equation system  */
 605 
 606       if(ido==-1){
 607         if(*isolver==0){
 608 #ifdef SPOOLES
 609           spooles_solve(temp_array,&neq[1]);
 610 #endif
 611         }
 612         else if(*isolver==4){
 613 #ifdef SGI
 614           sgi_solve(temp_array,token);
 615 #endif
 616         }
 617         else if(*isolver==5){
 618 #ifdef TAUCS
 619           tau_solve(temp_array,&neq[1]);
 620 #endif
 621         }
 622         else if(*isolver==7){
 623 #ifdef PARDISO
 624           pardiso_solve(temp_array,&neq[1],&symmetryflag);
 625 #endif
 626         }
 627         for(jrow=0;jrow<neq[1];jrow++){
 628           workd[ipntr[1]-1+jrow]=temp_array[jrow];
 629         }
 630       }
 631       else if(ido==1){
 632         if(*isolver==0){
 633 #ifdef SPOOLES
 634           spooles_solve(&workd[ipntr[2]-1],&neq[1]);
 635 #endif
 636         }
 637         else if(*isolver==4){
 638 #ifdef SGI
 639           sgi_solve(&workd[ipntr[2]-1],token);
 640 #endif
 641         }
 642         else if(*isolver==5){
 643 #ifdef TAUCS
 644           tau_solve(&workd[ipntr[2]-1],&neq[1]);
 645 #endif
 646         }
 647         else if(*isolver==7){
 648 #ifdef PARDISO
 649           pardiso_solve(&workd[ipntr[2]-1],&neq[1],&symmetryflag);
 650 #endif
 651         }
 652         for(jrow=0;jrow<neq[1];jrow++){
 653           workd[ipntr[1]-1+jrow]=workd[ipntr[2]-1+jrow];
 654         }
 655       }
 656 
 657     }
 658 
 659     if(ido==2){
 660         if(nasym==1){
 661             FORTRAN(opas,(&neq[1],&workd[ipntr[0]-1],&workd[ipntr[1]-1],
 662                     adb,aub,jq,irow,nzs));
 663         }else{
 664             FORTRAN(op,(&neq[1],&workd[ipntr[0]-1],&workd[ipntr[1]-1],
 665                     adb,aub,jq,irow));
 666         }
 667     }
 668 
 669     if(nasym==1){
 670         FORTRAN(dnaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,
 671                         iparam,ipntr,workd,workl,&lworkl,&info));
 672     }else{
 673         FORTRAN(dsaupd,(&ido,bmat,&neq[1],which,&nev,&tol,resid,&ncv,z,&ldz,
 674                 iparam,ipntr,workd,workl,&lworkl,&info));
 675     }
 676   }
 677 
 678 /*--------------------------------------------------------------------*/
 679 /*
 680    -----------
 681    free memory
 682    -----------
 683 */
 684   SFREE(temp_array);
 685   if(*isolver==0){
 686 #ifdef SPOOLES
 687     spooles_cleanup();
 688 #endif
 689   }
 690   else if(*isolver==4){
 691 #ifdef SGI
 692     sgi_cleanup(token);
 693 #endif
 694   }
 695   else if(*isolver==5){
 696 #ifdef TAUCS
 697     tau_cleanup();
 698 #endif
 699   }
 700   else if(*isolver==7){
 701 #ifdef PARDISO
 702       pardiso_cleanup(&neq[1],&symmetryflag);
 703 #endif
 704   }
 705 
 706   if(info!=0){
 707     printf("*ERROR in arpack: info=%" ITGFORMAT "\n",info);
 708     printf("       # of converged eigenvalues=%" ITGFORMAT "\n\n",iparam[4]);
 709   }         
 710 
 711   NNEW(select,ITG,ncv);
 712 
 713   if(nasym==1){
 714       NNEW(d,double,nev+1);
 715       NNEW(di,double,nev+1);
 716       NNEW(workev,double,3*ncv);
 717       FORTRAN(dneupd,(&rvec,howmny,select,d,di,z,&ldz,&sigma,&sigmai,
 718           workev,bmat,&neq[1],which,&nev,&tol,resid,
 719           &ncv,z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
 720       SFREE(workev);
 721       NNEW(dc,double,2*nev);
 722 
 723       /* storing as complex number and taking the square root */
 724 
 725       for(j=0;j<nev;j++){
 726           dc[2*j]=sqrt(sqrt(d[j]*d[j]+di[j]*di[j])+d[j])/sqrt(2.);
 727           dc[2*j+1]=sqrt(sqrt(d[j]*d[j]+di[j]*di[j])-d[j])/sqrt(2.);
 728           if(di[j]<0.) dc[2*j+1]=-dc[2*j+1];
 729       }
 730       FORTRAN(writeevcomplex,(dc,&nev,&fmin,&fmax));
 731       SFREE(di);SFREE(dc);
 732   }else{
 733       NNEW(d,double,nev);
 734       FORTRAN(dseupd,(&rvec,howmny,select,d,z,&ldz,&sigma,bmat,&neq[1],which,
 735           &nev,&tol,resid,&ncv,z,&ldz,iparam,ipntr,workd,workl,&lworkl,&info));
 736       FORTRAN(writeev,(d,&nev,&fmin,&fmax));
 737   }
 738   SFREE(select);SFREE(workd);SFREE(workl);SFREE(resid);
 739 
 740   /* writing the eigenvalues and mass matrix to a binary file */
 741 
 742   if(mei[3]==1){
 743 
 744       strcpy(fneig,jobnamec);
 745       strcat(fneig,".eig");
 746       
 747       if((f1=fopen(fneig,"wb"))==NULL){
 748           printf("*ERROR in arpack: cannot open eigenvalue file for writing...");
 749 
 750           exit(0);
 751       }
 752 
 753       /* storing a zero as indication that this was not a
 754          cyclic symmetry calculation */
 755 
 756       if(fwrite(&zero,sizeof(ITG),1,f1)!=1){
 757           printf("*ERROR saving the cyclic symmetry flag to the eigenvalue file...");
 758           exit(0);
 759       }
 760 
 761       /* Hermitian */
 762 
 763       if(fwrite(&nherm,sizeof(ITG),1,f1)!=1){
 764           printf("*ERROR saving the Hermitian flag to the eigenvalue file...");
 765           exit(0);
 766       }
 767 
 768       /* storing the number of eigenvalues */
 769 
 770       if(fwrite(&nev,sizeof(ITG),1,f1)!=1){
 771           printf("*ERROR saving the number of eigenvalues to the eigenvalue file...");
 772           exit(0);
 773       }
 774 
 775       /* the eigenfrequencies are stores as radians/time */
 776 
 777       if(fwrite(d,sizeof(double),nev,f1)!=nev){
 778           printf("*ERROR saving the eigenfrequencies to the eigenvalue file...");
 779           exit(0);
 780       }
 781 
 782       /* storing the stiffness matrix */
 783 
 784       if(fwrite(ad,sizeof(double),neq[1],f1)!=neq[1]){
 785           printf("*ERROR saving the diagonal of the stiffness matrix to the eigenvalue file...");
 786           exit(0);
 787       }
 788       if(fwrite(au,sizeof(double),nzs[2],f1)!=nzs[2]){
 789           printf("*ERROR saving the off-diagonal terms of the stiffness matrix to the eigenvalue file...");
 790           exit(0);
 791       }
 792 
 793       /* storing the mass matrix */
 794 
 795       if(fwrite(adb,sizeof(double),neq[1],f1)!=neq[1]){
 796           printf("*ERROR saving the diagonal of the mass matrix to the eigenvalue file...");
 797           exit(0);
 798       }
 799       if(fwrite(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
 800           printf("*ERROR saving the off-diagonal terms of the mass matrix to the eigenvalue file...");
 801           exit(0);
 802       }
 803   }
 804 
 805   /* calculating the participation factors and the relative effective
 806      modal mass */
 807 
 808   if(*ithermal!=2){
 809       FORTRAN(effectivemodalmass,(neq,nactdof,mi,adb,aub,jq,irow,&nev,z,co,nk));
 810   }
 811 
 812   SFREE(au);SFREE(ad);
 813 
 814   /* calculating the displacements and the stresses and storing */
 815   /* the results in frd format for each valid eigenmode */
 816 
 817   NNEW(v,double,mt**nk);
 818   NNEW(fn,double,mt**nk);
 819   NNEW(stn,double,6**nk);
 820   NNEW(inum,ITG,*nk);
 821   NNEW(stx,double,6*mi[0]**ne);
 822   if(*ithermal>1){
 823       NNEW(qfn,double,3**nk);
 824       NNEW(qfx,double,3*mi[0]**ne);
 825   }
 826   
 827   if(strcmp1(&filab[261],"E   ")==0) NNEW(een,double,6**nk);
 828   if(strcmp1(&filab[2697],"ME  ")==0) NNEW(emn,double,6**nk);
 829   if(strcmp1(&filab[522],"ENER")==0) NNEW(enern,double,*nk);
 830   if(strcmp1(&filab[2175],"CONT")==0) NNEW(cdn,double,6**nk);
 831 
 832   NNEW(temp_array,double,neq[1]);
 833 
 834   lfin=0;
 835   for(j=0;j<nev;++j){
 836     lint=lfin;
 837     lfin=lfin+neq[1];
 838 
 839     for(k=0;k<6*mi[0]*ne0;k++){eme[k]=0.;}
 840 
 841     sum=0.;
 842     for(k=0;k<neq[1];++k)
 843       temp_array[k]=0.;
 844     if(nasym==1){
 845         FORTRAN(opas,(&neq[1],&z[lint],temp_array,adb,aub,jq,irow,nzs));
 846     }else{
 847         FORTRAN(op,(&neq[1],&z[lint],temp_array,adb,aub,jq,irow));
 848     }
 849     for(k=0;k<neq[1];++k)
 850       sum+=z[lint+k]*temp_array[k];
 851     for(k=0;k<neq[1];++k)
 852       z[lint+k]=z[lint+k]/sqrt(sum);
 853 
 854     if(mei[3]==1){
 855         if(fwrite(&z[lint],sizeof(double),neq[1],f1)!=neq[1]){
 856             printf("*ERROR saving data to the eigenvalue file...");
 857             exit(0);
 858         }
 859     }
 860 
 861     /* check whether the frequency belongs to the requested
 862        interval */
 863 
 864     if(fmin>-0.5){
 865         if(fmin*fmin>d[j]) continue;
 866     }
 867     if(fmax>-0.5){
 868         if(fmax*fmax<d[j]) continue;
 869     }
 870 
 871     if(*nprint>0) FORTRAN(writehe,(&j));
 872 
 873     NNEW(eei,double,6*mi[0]*ne0);
 874     if(*nener==1){
 875         NNEW(stiini,double,6*mi[0]*ne0);
 876         NNEW(enerini,double,mi[0]*ne0);}
 877 
 878 //    memset(&v[0],0.,sizeof(double)*mt**nk);
 879     DMEMSET(v,0,mt**nk,0.);
 880     if(*iperturb==0){
 881       results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
 882               stx,elcon,
 883               nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
 884               norien,orab,ntmat_,t0,t0,ithermal,
 885               prestr,iprestr,filab,eme,emn,een,iperturb,
 886               f,fn,nactdof,&iout,qa,vold,&z[lint],
 887               nodeboun,ndirboun,xboun,nboun,ipompc,
 888               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
 889               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 890               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 891               &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
 892               emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
 893               iendset,
 894               ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 895               nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
 896               &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
 897               sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 898               mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 899               islavsurf,ielprop,prop);}
 900     else{
 901       results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
 902               stx,elcon,
 903               nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
 904               norien,orab,ntmat_,t0,t1old,ithermal,
 905               prestr,iprestr,filab,eme,emn,een,iperturb,
 906               f,fn,nactdof,&iout,qa,vold,&z[lint],
 907               nodeboun,ndirboun,xboun,nboun,ipompc,
 908               nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
 909               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 910               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
 911               &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
 912               xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
 913               ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 914               nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,&reltime,
 915               &ne0,xforc,nforc,thicke,shcon,nshcon,
 916               sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 917               mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 918               islavsurf,ielprop,prop);
 919     }
 920     SFREE(eei);
 921     if(*nener==1){
 922         SFREE(stiini);SFREE(enerini);}
 923 
 924     ++*kode;
 925     if(d[j]>=0.){
 926         freq=sqrt(d[j])/6.283185308;
 927     }else{
 928         freq=0.;
 929     }
 930 
 931     if(strcmp1(&filab[1044],"ZZS")==0){
 932         NNEW(neigh,ITG,40**ne);
 933         NNEW(ipneigh,ITG,*nk);
 934     }
 935 
 936     frd(co,nk,kon,ipkon,lakon,ne,v,stn,inum,nmethod,
 937             kode,filab,een,t1,fn,&freq,epn,ielmat,matname,enern,xstaten,
 938             nstate_,istep,&iinc,ithermal,qfn,&j,&noddiam,trab,inotr,
 939             ntrans,orab,ielorien,norien,description,ipneigh,neigh,
 940             mi,stx,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
 941             cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
 942             thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat);
 943 
 944     if(strcmp1(&filab[1044],"ZZS")==0){SFREE(ipneigh);SFREE(neigh);}
 945 
 946     if(*iperturb!=0){
 947         for(k=0;k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
 948             xstate[k]=xstateini[k];
 949         }         
 950     }
 951   }
 952 
 953   if((fmax>-0.5)&&(fmax*fmax>d[nev-1])){
 954     printf("\n*WARNING: not all frequencies in the requested interval might be found;\nincrease the number of requested frequencies\n");
 955   }
 956 
 957   if(mei[3]==1){
 958       fclose(f1);
 959   }
 960 
 961   if(*iperturb!=0){
 962       if(ncont!=0){
 963           *ne=ne0;*nkon=nkon0;
 964           if(*nener==1){
 965               RENEW(ener,double,mi[0]**ne*2);
 966           }
 967           RENEW(ipkon,ITG,*ne);
 968           RENEW(lakon,char,8**ne);
 969           RENEW(kon,ITG,*nkon);
 970           if(*norien>0){
 971               RENEW(ielorien,ITG,mi[2]**ne);
 972           }
 973           RENEW(ielmat,ITG,mi[2]**ne);
 974           SFREE(cg);SFREE(straight);
 975           SFREE(imastop);SFREE(itiefac);SFREE(islavnode);
 976           SFREE(nslavnode);SFREE(iponoels);SFREE(inoels);SFREE(imastnode);
 977           SFREE(nmastnode);SFREE(itietri);SFREE(koncont);SFREE(xnoels);
 978           SFREE(springarea);SFREE(xmastnor);
 979 
 980           if(*mortar==0){
 981               SFREE(areaslav);
 982           }else if(*mortar==1){
 983               SFREE(pmastsurf);SFREE(ipe);SFREE(ime);
 984               SFREE(islavact);
 985           }
 986       }
 987   }
 988   
 989   SFREE(adb);SFREE(aub);SFREE(temp_array);
 990 
 991   SFREE(v);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);
 992   SFREE(z);SFREE(d);SFREE(xstiff);SFREE(ipobody);
 993 
 994   if(*ithermal>1){SFREE(qfn);SFREE(qfx);}
 995 
 996   if(*nstate_!=0){SFREE(xstateini);}
 997 
 998   if(strcmp1(&filab[261],"E   ")==0) SFREE(een);
 999   if(strcmp1(&filab[2697],"ME  ")==0) SFREE(emn);
1000   if(strcmp1(&filab[522],"ENER")==0) SFREE(enern);
1001   if(strcmp1(&filab[2175],"CONT")==0) SFREE(cdn);
1002 
1003   for(k=0;k<6*mi[0]*ne0;k++){eme[k]=0.;}
1004 
1005   if(*iperturb!=0){
1006       mpcinfo[0]=memmpc_;mpcinfo[1]=mpcfree;mpcinfo[2]=icascade;
1007       mpcinfo[3]=maxlenmpc;
1008   }
1009 
1010   *irowp=irow;*enerp=ener;*xstatep=xstate;*ipkonp=ipkon;*lakonp=lakon;
1011   *konp=kon;*ielmatp=ielmat;*ielorienp=ielorien;
1012 
1013   *islavsurfp=islavsurf;*pslavsurfp=pslavsurf;*clearinip=clearini;
1014 
1015   return;
1016 }
1017 
1018 #endif

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