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 #include <unistd.h>
19 #include <stdio.h>
20 #include <math.h>
21 #include <stdlib.h>
22 #include <pthread.h>
23 #include "CalculiX.h"
24
25 static char *lakonf1;
26
27 static ITG num_cpus,*nef1,*ipnei1,*neifa1,*neiel1,*jq1,*irow1,*nzs1,
28 *ielfa1,*ifabou1,*nbody1,*neq1,*nactdohinv1,*iau61,*iturbulent1,
29 *inlet1;
30
31 static double *au1,*ad1,*b1,*vfa1,*xxn1,*area1,*vel1,
32 *umfa1,*alet1,*ale1,*gradkfa1,*xxi1,*body1,*volume1,*dtimef1,*velo1,
33 *veloo1,*cvfa1,*hcfa1,*cvel1,*gradvel1,*xload1,*xrlfa1,
34 *xxj1,*a11,*a21,*a31,*flux1,*xxni1,*xxnj1,*f11,*of21,*yy1,*umel1,
35 *gradkel1,*gradoel1,*sc1;
36
mafillkcompmain(ITG * nef,ITG * ipnei,ITG * neifa,ITG * neiel,double * vfa,double * xxn,double * area,double * au,double * ad,ITG * jq,ITG * irow,ITG * nzs,double * b,double * vel,double * umfa,double * alet,double * ale,double * gradkfa,double * xxi,double * body,double * volume,ITG * ielfa,char * lakonf,ITG * ifabou,ITG * nbody,ITG * neq,double * dtimef,double * velo,double * veloo,double * cvfa,double * hcfa,double * cvel,double * gradvel,double * xload,double * xrlfa,double * xxj,ITG * nactdohinv,double * a1,double * a2,double * a3,double * flux,ITG * iau6,double * xxni,double * xxnj,ITG * iturbulent,double * f1,double * of2,double * yy,double * umel,double * gradkel,double * gradoel,ITG * inlet,double * sc)37 void mafillkcompmain(ITG *nef,ITG *ipnei,ITG *neifa,
38 ITG *neiel,double *vfa,double *xxn,double *area,
39 double *au,double *ad,ITG *jq,ITG *irow,ITG *nzs,
40 double *b,double *vel,double *umfa,double *alet,
41 double *ale,double *gradkfa,double *xxi,double *body,
42 double *volume,ITG *ielfa,char *lakonf,
43 ITG *ifabou,ITG *nbody,ITG *neq,double *dtimef,double *velo,
44 double *veloo,double *cvfa,double *hcfa,double *cvel,
45 double *gradvel,double *xload,double *xrlfa,
46 double *xxj,ITG *nactdohinv,double *a1,double *a2,double *a3,
47 double *flux,ITG *iau6,double *xxni,double *xxnj,
48 ITG *iturbulent,double *f1,double *of2,double *yy,
49 double *umel,double *gradkel,double *gradoel,
50 ITG *inlet,double *sc){
51
52 ITG i;
53
54 /* variables for multithreading procedure */
55
56 ITG sys_cpus,*ithread=NULL;
57 char *env,*envloc,*envsys;
58
59 num_cpus = 0;
60 sys_cpus=0;
61
62 /* explicit user declaration prevails */
63
64 envsys=getenv("NUMBER_OF_CPUS");
65 if(envsys){
66 sys_cpus=atoi(envsys);
67 if(sys_cpus<0) sys_cpus=0;
68 }
69
70 /* automatic detection of available number of processors */
71
72 if(sys_cpus==0){
73 sys_cpus = getSystemCPUs();
74 if(sys_cpus<1) sys_cpus=1;
75 }
76
77 /* local declaration prevails, if strictly positive */
78
79 envloc = getenv("CCX_NPROC_CFD");
80 if(envloc){
81 num_cpus=atoi(envloc);
82 if(num_cpus<0){
83 num_cpus=0;
84 }else if(num_cpus>sys_cpus){
85 num_cpus=sys_cpus;
86 }
87
88 }
89
90 /* else global declaration, if any, applies */
91
92 env = getenv("OMP_NUM_THREADS");
93 if(num_cpus==0){
94 if (env)
95 num_cpus = atoi(env);
96 if (num_cpus < 1) {
97 num_cpus=1;
98 }else if(num_cpus>sys_cpus){
99 num_cpus=sys_cpus;
100 }
101 }
102
103 // next line is to be inserted in a similar way for all other paralell parts
104
105 if(*nef<num_cpus) num_cpus=*nef;
106
107 pthread_t tid[num_cpus];
108
109 /* calculating the stiffness and/or mass matrix
110 (symmetric part) */
111
112 nef1=nef;ipnei1=ipnei;neifa1=neifa;neiel1=neiel;vfa1=vfa;xxn1=xxn;
113 area1=area;
114 jq1=jq;irow1=irow;nzs1=nzs;vel1=vel;umfa1=umfa;alet1=alet;ale1=ale;
115 gradkfa1=gradkfa;xxi1=xxi;body1=body;volume1=volume;
116 ielfa1=ielfa;lakonf1=lakonf;ifabou1=ifabou;
117 nbody1=nbody;neq1=neq;dtimef1=dtimef;velo1=velo;veloo1=veloo;
118 cvfa1=cvfa;hcfa1=hcfa;cvel1=cvel;gradvel1=gradvel;xload1=xload;
119 xrlfa1=xrlfa;xxj1=xxj;nactdohinv1=nactdohinv;a11=a1;
120 a21=a2;a31=a3;flux1=flux;iau61=iau6;ad1=ad;au1=au;b1=b;xxni1=xxni;
121 xxnj1=xxnj,iturbulent1=iturbulent,f11=f1,of21=of2;yy1=yy;umel1=umel;
122 gradkel1=gradkel;gradoel1=gradoel;inlet1=inlet;sc1=sc;
123
124 /* create threads and wait */
125
126 NNEW(ithread,ITG,num_cpus);
127 for(i=0; i<num_cpus; i++) {
128 ithread[i]=i;
129 pthread_create(&tid[i], NULL, (void *)mafillkcompmt, (void *)&ithread[i]);
130 }
131 for(i=0; i<num_cpus; i++) pthread_join(tid[i], NULL);
132
133 SFREE(ithread);
134
135 return;
136
137 }
138
139 /* subroutine for multithreading of mafillkcomp */
140
mafillkcompmt(ITG * i)141 void *mafillkcompmt(ITG *i){
142
143 ITG nefa,nefb,nefdelta;
144
145 // ceil -> floor
146
147 nefdelta=(ITG)floor(*nef1/(double)num_cpus);
148 nefa=*i*nefdelta+1;
149 nefb=(*i+1)*nefdelta;
150 // next line! -> all parallel sections
151 if((*i==num_cpus-1)&&(nefb<*nef1)) nefb=*nef1;
152
153 FORTRAN(mafillkcomp,(nef1,ipnei1,neifa1,neiel1,vfa1,xxn1,area1,
154 au1,ad1,jq1,irow1,nzs1,
155 b1,vel1,umfa1,alet1,ale1,gradkfa1,xxi1,
156 body1,volume1,ielfa1,lakonf1,ifabou1,
157 nbody1,neq1,dtimef1,velo1,veloo1,cvfa1,hcfa1,cvel1,
158 gradvel1,xload1,xrlfa1,xxj1,nactdohinv1,
159 a11,a21,a31,flux1,&nefa,&nefb,iau61,xxni1,xxnj1,
160 iturbulent1,f11,of21,yy1,umel1,gradkel1,gradoel1,
161 inlet1,sc1));
162
163 return NULL;
164 }
165