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 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 #ifdef PASTIX
41 #include "pastix.h"
42 #endif
43
44
dudsmain(ITG * isolver,double * au,double * ad,double * aub,double * adb,ITG * icol,ITG * irow,ITG * jq,ITG * neq,ITG * nzs,double * df,ITG * jqs,ITG * irows,ITG * ndesi,double * duds)45 void dudsmain(ITG *isolver,double *au,double *ad,double *aub,double*adb,
46 ITG *icol,ITG *irow,ITG *jq,ITG *neq,ITG *nzs,double *df,ITG *jqs,
47 ITG *irows,ITG *ndesi,double *duds){
48
49 /* calculating the matrix duds */
50
51 ITG symmetryflag=0,inputformat=0,idesvar,idof,j,nrhs=1,num_cpus;
52
53 double sigma=0,*dudsvec=NULL;
54
55 /* variables for multithreading procedure */
56
57 ITG sys_cpus,*ithread=NULL;
58 char *env,*envloc,*envsys;
59
60 num_cpus=0;
61 sys_cpus=0;
62
63 /* explicit user declaration prevails */
64
65 envsys=getenv("NUMBER_OF_CPUS");
66 if(envsys){
67 sys_cpus=atoi(envsys);
68 if(sys_cpus<0) sys_cpus=0;
69 }
70
71 /* automatic detection of available number of processors */
72
73 if(sys_cpus==0){
74 sys_cpus = getSystemCPUs();
75 if(sys_cpus<1) sys_cpus=1;
76 }
77
78 /* local declaration prevails, if strictly positive */
79
80 envloc = getenv("CCX_NPROC_SENS");
81 if(envloc){
82 num_cpus=atoi(envloc);
83 if(num_cpus<0){
84 num_cpus=0;
85 }else if(num_cpus>sys_cpus){
86 num_cpus=sys_cpus;
87 }
88
89 }
90
91 /* else global declaration, if any, applies */
92
93 env = getenv("OMP_NUM_THREADS");
94 if(num_cpus==0){
95 if (env)
96 num_cpus = atoi(env);
97 if (num_cpus < 1) {
98 num_cpus=1;
99 }else if(num_cpus>sys_cpus){
100 num_cpus=sys_cpus;
101 }
102 }
103
104 pthread_t tid[num_cpus];
105
106 /* factorize the system */
107
108 if(*isolver==0){
109 #ifdef SPOOLES
110 spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
111 &symmetryflag,&inputformat,&nzs[2]);
112 #else
113 printf("*ERROR in dudsmain: the SPOOLES library is not linked\n\n");
114 FORTRAN(stop,());
115 #endif
116 }
117 else if(*isolver==4){
118 #ifdef SGI
119 token=1;
120 sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
121 #else
122 printf("*ERROR in dudsmain: the SGI library is not linked\n\n");
123 FORTRAN(stop,());
124 #endif
125 }
126 else if(*isolver==5){
127 #ifdef TAUCS
128 tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
129 #else
130 printf("*ERROR in dudsmain: the TAUCS library is not linked\n\n");
131 FORTRAN(stop,());
132 #endif
133 }
134 else if(*isolver==7){
135 #ifdef PARDISO
136 pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
137 &symmetryflag,&inputformat,jq,&nzs[2]);
138 #else
139 printf("*ERROR in dudsmain: the PARDISO library is not linked\n\n");
140 FORTRAN(stop,());
141 #endif
142 }
143 else if(*isolver==8){
144 #ifdef PASTIX
145 pastix_factor_main(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
146 &symmetryflag,&inputformat,jq,&nzs[2]);
147 #else
148 printf("*ERROR in dudsmain: the PASTIX library is not linked\n\n");
149 FORTRAN(stop,());
150 #endif
151 }
152
153 /* Computation of the matrix duds */
154 /* duds = K^-1 * ( dF/ds + dK/ds * u )
155 /* = K^-1 * df */
156 /* ndesi system of equations have to be solved to determine duds */
157
158 NNEW(dudsvec,double,neq[1]);
159
160 for(idesvar=0;idesvar<*ndesi;idesvar++){
161
162 /* initialize vector dudsvec */
163
164 for(j=0;j<neq[1];j++){
165
166 dudsvec[j]=0.;
167
168 }
169
170 /* assembly of the vector taken from df */
171
172 for(j=jqs[idesvar];j<jqs[idesvar+1]-1;j++){
173
174 idof=irows[j-1]-1;
175 dudsvec[idof]=df[j-1];
176
177 }
178
179 /*Alternative way to assemble the vector dudsvec /*
180 /*FORTRAN(dudsassembly,(jqs,irows,neq,df,dudsvec,&idesvar));*/
181
182 /* solve the system */
183
184 if(*isolver==0){
185 #ifdef SPOOLES
186 spooles_solve(dudsvec,&neq[1]);
187 #endif
188 }
189 else if(*isolver==4){
190 #ifdef SGI
191 sgi_solve(dudsvec,token);
192 #endif
193 }
194 else if(*isolver==5){
195 #ifdef TAUCS
196 tau_solve(dudsvec,&neq[1]);
197 #endif
198 }
199 else if(*isolver==7){
200 #ifdef PARDISO
201 pardiso_solve(dudsvec,&neq[1],&symmetryflag,&inputformat,&nrhs);
202 #endif
203 }
204 else if(*isolver==8){
205 #ifdef PASTIX
206 pastix_solve(dudsvec,&neq[1],&symmetryflag,&nrhs);
207 #endif
208 }
209
210 /* copy results vector in duds */
211
212 for(j=0;j<neq[1];j++){
213
214 duds[(neq[1]*idesvar)+j]=dudsvec[j];
215
216 }
217
218 }
219
220 SFREE(dudsvec);
221
222 /* clean the system */
223
224 if(*isolver==0){
225 #ifdef SPOOLES
226 spooles_cleanup();
227 #endif
228 }
229 else if(*isolver==4){
230 #ifdef SGI
231 sgi_cleanup(token);
232 #endif
233 }
234 else if(*isolver==5){
235 #ifdef TAUCS
236 tau_cleanup();
237 #endif
238 }
239 else if(*isolver==7){
240 #ifdef PARDISO
241 pardiso_cleanup(&neq[1],&symmetryflag,&inputformat);
242 #endif
243 }
244 else if(*isolver==8){
245 #ifdef PASTIX
246 #endif
247 }
248
249 return;
250
251 }
252
253 #endif
254