root/src/pardiso_as.c

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

DEFINITIONS

This source file includes following definitions.
  1. pardiso_factor_as
  2. pardiso_solve_as
  3. pardiso_cleanup_as
  4. pardiso_main_as

   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 *irowpardisoas=NULL,*pointersas=NULL,iparmas[64];
  27 long long ptas[64];
  28 double *aupardisoas=NULL;
  29 ITG nthread_mkl_as=0;
  30 
  31 /* factorization, solving and cleaning with PARDISO for
  32    real unsymmetric matrices */
  33 
  34 void pardiso_factor_as(double *ad, double *au, double *adb, double *aub, 
  35                 double *sigma,ITG *icol, ITG *irow, 
  36                 ITG *neq, ITG *nzs, ITG *jq){
  37 
  38   char *env;
  39   ITG i,j,k,maxfct=1,mnum=1,mtype=11,phase=12,nrhs=1,*perm=NULL,
  40       msglvl=0,error=0,ifortran,lfortran,index,id;
  41   long long ndim;
  42   ITG nthread,nthread_v;
  43   double *b=NULL,*x=NULL;
  44 
  45   printf(" Factoring the system of equations using the asymmetric pardiso solver\n");
  46 
  47   iparmas[0]=0;
  48 /* set MKL_NUM_THREADS to min(CCX_NPROC_EQUATION_SOLVER,OMP_NUM_THREADS)
  49    must be done once  */
  50   if (nthread_mkl_as == 0) {
  51     nthread=1;
  52     env=getenv("MKL_NUM_THREADS");
  53     if(env) {
  54       nthread=atoi(env);}
  55     else {
  56       env=getenv("OMP_NUM_THREADS");
  57       if(env) {nthread=atoi(env);}
  58     }
  59     env=getenv("CCX_NPROC_EQUATION_SOLVER");
  60     if(env) {
  61       nthread_v=atoi(env);
  62       if (nthread_v <= nthread) {nthread=nthread_v;}
  63     }
  64     if (nthread < 1) {nthread=1;}
  65     sprintf(envMKL,"MKL_NUM_THREADS=%" ITGFORMAT "",nthread);  
  66     putenv(envMKL);
  67     nthread_mkl_as=nthread;
  68   }
  69     
  70   printf(" number of threads =% d\n\n",nthread_mkl_as);
  71 
  72   for(i=0;i<64;i++){ptas[i]=0;}
  73 
  74   ndim=*neq+(long long)2**nzs;
  75 
  76   NNEW(pointersas,ITG,*neq+1);
  77   NNEW(irowpardisoas,ITG,ndim);
  78   NNEW(aupardisoas,double,ndim);
  79 
  80   k=0;
  81 
  82   /* PARDISO requires the matrix to be stored row by row
  83      aupardisoas contains the entries
  84      irowpardisoas the corresponding column numbers
  85          (for each row in ascending order)
  86      pointersas(i) points to the first entry for row i in 
  87            field aupardisoas */
  88 
  89   if(*sigma==0.){
  90     for(i=0;i<*neq;i++){
  91       pointersas[i]=k+1;
  92       
  93       /* lower left triangular matrix */
  94 
  95       for(j=0;j<i;j++){
  96         ifortran=i+1;
  97         lfortran=jq[j+1]-jq[j];
  98         FORTRAN(nident,(&irow[jq[j]-1],&ifortran,&lfortran,&id));
  99         if(id>0){
 100           index=jq[j]+id-2;
 101           if(irow[index]==ifortran){
 102             irowpardisoas[k]=j+1;
 103             aupardisoas[k]=au[index];
 104             k++;
 105           }
 106         }
 107       }
 108 
 109       /* diagonal entry */
 110 
 111       irowpardisoas[k]=i+1;
 112       aupardisoas[k]=ad[i];
 113       k++;
 114 
 115       /* upper right triangular matrix */
 116 
 117       for(j=jq[i];j<jq[i+1];j++){
 118         irowpardisoas[k]=irow[j-1];
 119         aupardisoas[k]=au[j+*nzs-1];
 120         k++;
 121       }
 122     }
 123     pointersas[*neq]=k+1;
 124   }else{
 125     for(i=0;i<*neq;i++){
 126       pointersas[i]=k+1;
 127       
 128       /* lower left triangular matrix */
 129 
 130       for(j=0;j<i;j++){
 131         ifortran=i+1;
 132         lfortran=jq[j+1]-jq[j];
 133         FORTRAN(nident,(&irow[jq[j]-1],&ifortran,&lfortran,&id));
 134         if(id>0){
 135           index=jq[j]+id-2;
 136           if(irow[index]==ifortran){
 137             irowpardisoas[k]=j+1;
 138             aupardisoas[k]=au[index]-*sigma*aub[index];
 139             k++;
 140           }
 141         }
 142       }
 143 
 144       /* diagonal entry */
 145 
 146       irowpardisoas[k]=i+1;
 147       aupardisoas[k]=ad[i]-*sigma*adb[i];
 148       k++;
 149 
 150       /* upper right triangular matrix */
 151 
 152       for(j=jq[i];j<jq[i+1];j++){
 153         irowpardisoas[k]=irow[j-1];
 154         aupardisoas[k]=au[j+*nzs-1]-*sigma*aub[j+*nzs-1];
 155         k++;
 156       }
 157     }
 158     pointersas[*neq]=k+1;
 159   }
 160 
 161   FORTRAN(pardiso,(ptas,&maxfct,&mnum,&mtype,&phase,neq,aupardisoas,
 162                    pointersas,irowpardisoas,perm,&nrhs,iparmas,&msglvl,
 163                    b,x,&error));
 164 
 165   return;
 166 }
 167 
 168 void pardiso_solve_as(double *b, ITG *neq){
 169 
 170   char *env;
 171   ITG maxfct=1,mnum=1,mtype=11,phase=33,*perm=NULL,nrhs=1,
 172       msglvl=0,i,error=0;
 173   double *x=NULL;
 174 
 175   printf(" Solving the system of equations using the asymmetric pardiso solver\n");
 176 
 177   iparmas[0]=0;
 178 /* pardiso_factor has been called befor, MKL_NUM_THREADS=nthread_mkl_as is set*/
 179 
 180   printf(" number of threads =% d\n\n",nthread_mkl_as);
 181 
 182   NNEW(x,double,*neq);
 183 
 184   FORTRAN(pardiso,(ptas,&maxfct,&mnum,&mtype,&phase,neq,aupardisoas,
 185                    pointersas,irowpardisoas,perm,&nrhs,iparmas,&msglvl,
 186                    b,x,&error));
 187 
 188   for(i=0;i<*neq;i++){b[i]=x[i];}
 189   SFREE(x);
 190 
 191   return;
 192 }
 193 
 194 void pardiso_cleanup_as(ITG *neq){
 195 
 196   ITG maxfct=1,mnum=1,mtype=11,phase=-1,*perm=NULL,nrhs=1,
 197       msglvl=0,error=0;
 198   double *b=NULL,*x=NULL;
 199 
 200   FORTRAN(pardiso,(ptas,&maxfct,&mnum,&mtype,&phase,neq,aupardisoas,
 201                    pointersas,irowpardisoas,perm,&nrhs,iparmas,&msglvl,
 202                    b,x,&error));
 203 
 204   SFREE(irowpardisoas);
 205   SFREE(aupardisoas);
 206   SFREE(pointersas);
 207 
 208   return;
 209 }
 210 
 211 void pardiso_main_as(double *ad, double *au, double *adb, double *aub, double *sigma,
 212          double *b, ITG *icol, ITG *irow, 
 213          ITG *neq, ITG *nzs, ITG *jq){
 214 
 215   if(*neq==0) return;
 216 
 217   pardiso_factor_as(ad,au,adb,aub,sigma,icol,irow, 
 218              neq,nzs,jq);
 219 
 220   pardiso_solve_as(b,neq);
 221 
 222   pardiso_cleanup_as(neq);
 223 
 224   return;
 225 }
 226 
 227 #endif
 228 

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