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