root/src/pardiso.c

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

DEFINITIONS

This source file includes following definitions.
  1. pardiso_factor
  2. pardiso_solve
  3. pardiso_cleanup
  4. pardiso_main

   1 /*     CalculiX - A 3-dimensional finite element program                   */
   2 /*              Copyright (C) 1998-2005 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 PARDISO
  19 
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23 #include "CalculiX.h"
  24 #include "pardiso.h"
  25 
  26 ITG *icolpardiso=NULL,*pointers=NULL,iparm[64];
  27 long long pt[64];
  28 double *aupardiso=NULL;
  29 /* double dparm[64];  not used */
  30 ITG nthread_mkl=0;
  31 /* char envMKL[32];   moved to pardiso.h */
  32 
  33 void pardiso_factor(double *ad, double *au, double *adb, double *aub, 
  34                 double *sigma,ITG *icol, ITG *irow, 
  35                 ITG *neq, ITG *nzs, ITG *symmetryflag, ITG *inputformat,
  36                 ITG *jq, ITG *nzs3){
  37 
  38   char *env;
  39 /*  char env1[32]; */
  40   ITG i,j,k,l,maxfct=1,mnum=1,phase=12,nrhs=1,*perm=NULL,mtype,
  41       msglvl=0,error=0,*irowpardiso=NULL,kflag,kstart,n,ifortran,
  42       lfortran,index,id,k2;
  43   ITG ndim,nthread,nthread_v;
  44   double *b=NULL,*x=NULL;
  45 
  46   if(*symmetryflag==0){
  47       printf(" Factoring the system of equations using the symmetric pardiso solver\n");
  48   }else{
  49       printf(" Factoring the system of equations using the unsymmetric pardiso solver\n");
  50   }
  51 
  52   iparm[0]=0;
  53 /* set MKL_NUM_THREADS to min(CCX_NPROC_EQUATION_SOLVER,OMP_NUM_THREADS)
  54    must be done once  */
  55   if (nthread_mkl == 0) {
  56     nthread=1;
  57     env=getenv("MKL_NUM_THREADS");
  58     if(env) {
  59       nthread=atoi(env);}
  60     else {
  61       env=getenv("OMP_NUM_THREADS");
  62       if(env) {nthread=atoi(env);}
  63     }
  64     env=getenv("CCX_NPROC_EQUATION_SOLVER");
  65     if(env) {
  66       nthread_v=atoi(env);
  67       if (nthread_v <= nthread) {nthread=nthread_v;}
  68     }
  69     if (nthread < 1) {nthread=1;}
  70     sprintf(envMKL,"MKL_NUM_THREADS=%" ITGFORMAT "",nthread);  
  71     putenv(envMKL);
  72     nthread_mkl=nthread;
  73   }
  74     
  75   printf(" number of threads =% d\n\n",nthread_mkl);
  76 
  77   for(i=0;i<64;i++){pt[i]=0;}
  78 
  79   if(*symmetryflag==0){
  80 
  81       /* symmetric matrix; the subdiagonal entries are stored
  82          column by column in au, the diagonal entries in ad;
  83          pardiso needs the entries row per row */      
  84 
  85       mtype=-2;
  86       
  87       ndim=*neq+*nzs;
  88       
  89       NNEW(pointers,ITG,*neq+1);
  90       NNEW(icolpardiso,ITG,ndim);
  91       NNEW(aupardiso,double,ndim);
  92       
  93       k=ndim;
  94       l=*nzs;
  95       
  96       if(*sigma==0.){
  97           pointers[*neq]=ndim+1;
  98           for(i=*neq-1;i>=0;--i){
  99               for(j=0;j<icol[i];++j){
 100                   icolpardiso[--k]=irow[--l];
 101                   aupardiso[k]=au[l];
 102               }
 103               pointers[i]=k--;
 104               icolpardiso[k]=i+1;
 105               aupardiso[k]=ad[i];
 106           }
 107       }
 108       else{
 109           pointers[*neq]=ndim+1;
 110           for(i=*neq-1;i>=0;--i){
 111               for(j=0;j<icol[i];++j){
 112                   icolpardiso[--k]=irow[--l];
 113                   aupardiso[k]=au[l]-*sigma*aub[l];
 114               }
 115               pointers[i]=k--;
 116               icolpardiso[k]=i+1;
 117               aupardiso[k]=ad[i]-*sigma*adb[i];
 118           }
 119       }
 120   }else{
 121       mtype=11;
 122 
 123       if(*inputformat==3){
 124 
 125           /* off-diagonal terms  are stored column per
 126              column from top to bottom in au;
 127              diagonal terms are stored in ad  */
 128 
 129           ndim=*neq+*nzs;
 130           NNEW(pointers,ITG,*neq+1);
 131           NNEW(irowpardiso,ITG,ndim);     
 132           NNEW(icolpardiso,ITG,ndim);
 133           NNEW(aupardiso,double,ndim);
 134           
 135           k=0;
 136           k2=0;
 137           for(i=0;i<*neq;i++){
 138               for(j=0;j<icol[i];j++){
 139                   if(au[k]>1.e-12||au[k]<-1.e-12){
 140                       icolpardiso[k2]=i+1;
 141                       irowpardiso[k2]=irow[k];
 142                       aupardiso[k2]=au[k];
 143                       k2++;               
 144                   }
 145                   k++;        
 146               }   
 147           }  
 148           /* diagonal terms */  
 149           for(i=0;i<*neq;i++){
 150               icolpardiso[k2]=i+1;
 151               irowpardiso[k2]=i+1;
 152               aupardiso[k2]=ad[i];
 153               k2++;       
 154           }
 155           ndim=k2;
 156           
 157           /* pardiso needs the entries row per row; so sorting is
 158              needed */ 
 159           
 160           kflag=2;
 161           FORTRAN(isortiid,(irowpardiso,icolpardiso,aupardiso,
 162                             &ndim,&kflag));
 163           
 164           /* sorting each row */
 165           
 166           k=0;
 167           pointers[0]=1;
 168           for(i=0;i<*neq;i++){
 169               j=i+1;
 170               kstart=k;
 171               do{
 172                   if(irowpardiso[k]!=j ){
 173                       n=k-kstart;                 
 174                       FORTRAN(isortid,(&icolpardiso[kstart],&aupardiso[kstart],
 175                                        &n,&kflag));
 176                       pointers[i+1]=k+1;
 177                       break;  
 178                   }else{
 179                       if(k+1==ndim){
 180                           n=k-kstart+1;   
 181                           FORTRAN(isortid,(&icolpardiso[kstart],
 182                                   &aupardiso[kstart],&n,&kflag));
 183                           break;               
 184                       }else{
 185                           k++;         
 186                       }  
 187                   }
 188               }while(1);
 189           }
 190           pointers[*neq]=ndim+1;
 191           SFREE(irowpardiso);
 192 
 193       }else if(*inputformat==1){
 194           
 195           /* lower triangular matrix is stored column by column in
 196              au, followed by the upper triangular matrix row by row;
 197              the diagonal terms are stored in ad */
 198 
 199           /* reordering lower triangular matrix */
 200 
 201           ndim=*nzs;
 202           NNEW(pointers,ITG,*neq+1);
 203           NNEW(irowpardiso,ITG,ndim);
 204           NNEW(icolpardiso,ITG,ndim);
 205           NNEW(aupardiso,double,ndim);
 206           
 207           k=0;
 208           for(i=0;i<*neq;i++){
 209               for(j=0;j<icol[i];j++){
 210                   icolpardiso[k]=i+1;
 211                   irowpardiso[k]=irow[k];
 212                   aupardiso[k]=au[k];
 213                   k++;
 214               }
 215           }
 216           
 217           /* pardiso needs the entries row per row; so sorting is
 218              needed */
 219           
 220           kflag=2;
 221           FORTRAN(isortiid,(irowpardiso,icolpardiso,aupardiso,
 222           &ndim,&kflag));
 223           
 224           /* sorting each row */
 225           
 226           k=0;
 227           pointers[0]=1;
 228           for(i=0;i<*neq;i++){
 229               j=i+1;
 230               kstart=k;
 231               do{
 232 
 233                   /* end of row reached */
 234 
 235                   if(irowpardiso[k]!=j){
 236                       n=k-kstart;
 237                       FORTRAN(isortid,(&icolpardiso[kstart],&aupardiso[kstart],
 238                                        &n,&kflag));
 239                       pointers[i+1]=k+1;
 240                       break;
 241                   }else{
 242 
 243                       /* end of last row reached */
 244 
 245                       if(k+1==ndim){
 246                           n=k-kstart+1;
 247                           FORTRAN(isortid,(&icolpardiso[kstart],&aupardiso[kstart],
 248                                            &n,&kflag));
 249                           break;
 250                       }else{
 251 
 252                           /* end of row not yet reached */
 253 
 254                           k++;
 255                       }
 256                   }
 257               }while(1);
 258           }
 259           pointers[*neq]=ndim+1;
 260           SFREE(irowpardiso);
 261 
 262           /* composing the matrix: lower triangle + diagonal + upper triangle */
 263 
 264           ndim=*neq+2**nzs;
 265           RENEW(icolpardiso,ITG,ndim);
 266           RENEW(aupardiso,double,ndim);
 267           k=ndim;
 268           for(i=*neq-1;i>=0;i--){
 269               l=k+1;
 270               for(j=jq[i+1]-1;j>=jq[i];j--){
 271                   icolpardiso[--k]=irow[j-1];
 272                   aupardiso[k]=au[j+*nzs3-1];
 273               }
 274               icolpardiso[--k]=i+1;
 275               aupardiso[k]=ad[i];
 276               for(j=pointers[i+1]-1;j>=pointers[i];j--){
 277                   icolpardiso[--k]=icolpardiso[j-1];
 278                   aupardiso[k]=aupardiso[j-1];
 279               }
 280               pointers[i+1]=l;
 281           }
 282           pointers[0]=1;
 283       }
 284   }
 285 
 286   FORTRAN(pardiso,(pt,&maxfct,&mnum,&mtype,&phase,neq,aupardiso,
 287                    pointers,icolpardiso,perm,&nrhs,iparm,&msglvl,
 288                    b,x,&error));
 289 
 290   return;
 291 }
 292 
 293 void pardiso_solve(double *b, ITG *neq,ITG *symmetryflag){
 294 
 295 /*  char *env;
 296   char env1[32]; 
 297   ITG nthread,nthread_v; */
 298   ITG maxfct=1,mnum=1,phase=33,*perm=NULL,nrhs=1,mtype,
 299       msglvl=0,i,error=0;
 300   double *x=NULL;
 301 
 302   if(*symmetryflag==0){
 303       printf(" Solving the system of equations using the symmetric pardiso solver\n");
 304   }else{
 305       printf(" Solving the system of equations using the unsymmetric pardiso solver\n");
 306   }
 307 
 308   if(*symmetryflag==0){
 309       mtype=-2;
 310   }else{
 311       mtype=11;
 312   }
 313   iparm[0]=0;
 314 /* pardiso_factor has been called befor, MKL_NUM_THREADS=nthread_mkl is set*/
 315 
 316   printf(" number of threads =% d\n\n",nthread_mkl);
 317 
 318   NNEW(x,double,*neq);
 319 
 320   FORTRAN(pardiso,(pt,&maxfct,&mnum,&mtype,&phase,neq,aupardiso,
 321                    pointers,icolpardiso,perm,&nrhs,iparm,&msglvl,
 322                    b,x,&error));
 323 
 324   for(i=0;i<*neq;i++){b[i]=x[i];}
 325   SFREE(x);
 326 
 327   return;
 328 }
 329 
 330 void pardiso_cleanup(ITG *neq,ITG *symmetryflag){
 331 
 332   ITG maxfct=1,mnum=1,phase=-1,*perm=NULL,nrhs=1,mtype,
 333       msglvl=0,error=0;
 334   double *b=NULL,*x=NULL;
 335 
 336   if(*symmetryflag==0){
 337       mtype=-2;
 338   }else{
 339       mtype=11;
 340   }
 341 
 342   FORTRAN(pardiso,(pt,&maxfct,&mnum,&mtype,&phase,neq,aupardiso,
 343                    pointers,icolpardiso,perm,&nrhs,iparm,&msglvl,
 344                    b,x,&error));
 345 
 346   SFREE(icolpardiso);
 347   SFREE(aupardiso);
 348   SFREE(pointers);
 349 
 350   return;
 351 }
 352 
 353 void pardiso_main(double *ad, double *au, double *adb, double *aub, 
 354          double *sigma,double *b, ITG *icol, ITG *irow, 
 355          ITG *neq, ITG *nzs,ITG *symmetryflag,ITG *inputformat,
 356          ITG *jq, ITG *nzs3){
 357 
 358   if(*neq==0) return;
 359 
 360   pardiso_factor(ad,au,adb,aub,sigma,icol,irow, 
 361                  neq,nzs,symmetryflag,inputformat,jq,nzs3);
 362 
 363   pardiso_solve(b,neq,symmetryflag);
 364 
 365   pardiso_cleanup(neq,symmetryflag);
 366 
 367   return;
 368 }
 369 
 370 #endif
 371 

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