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
randomfieldmain(ITG * kon,ITG * ipkon,char * lakon,ITG * ne,ITG * nmpc,ITG * nactdof,ITG * mi,ITG * nodedesi,ITG * ndesi,ITG * istartdesi,ITG * ialdesi,double * co,double * physcon,ITG * isolver,ITG * ntrans,ITG * nk,ITG * inotr,double * trab,char * jobnamec,ITG * nboun,double * cs,ITG * mcs,ITG * inum,ITG * nmethod,ITG * kode,char * filab,ITG * nstate_,ITG * istep,char * description,char * set,ITG * nset,ITG * iendset,char * output,ITG * istartset,ITG * ialset,double * extnor,ITG * irandomtype,double * randomval,ITG * irobustdesign,ITG * ndesibou,ITG * nodedesibou,ITG * nodedesiinvbou)45 void randomfieldmain(ITG *kon,ITG *ipkon,char *lakon,ITG *ne,ITG *nmpc,
46 ITG *nactdof,ITG *mi,ITG *nodedesi,ITG *ndesi,
47 ITG *istartdesi,
48 ITG *ialdesi,double *co,double *physcon,ITG *isolver,
49 ITG *ntrans,
50 ITG *nk,ITG *inotr,double *trab,char *jobnamec,ITG *nboun,
51 double *cs,
52 ITG *mcs,ITG *inum,ITG *nmethod,ITG *kode,char *filab,
53 ITG *nstate_,
54 ITG *istep,char *description,char *set,ITG *nset,
55 ITG *iendset,
56 char *output,ITG *istartset,ITG *ialset,double *extnor,
57 ITG *irandomtype,double *randomval,ITG *irobustdesign,
58 ITG *ndesibou,
59 ITG *nodedesibou,ITG *nodedesiinvbou){
60
61 /* generating the random field for robust assessment */
62
63 char bmat[2]="I", which[3]="SA", howmny[2]="A",*objectset=NULL,*orname=NULL;
64
65 ITG nzss,nzsd,*mast1=NULL,*irows=NULL,*icols=NULL,*jqs=NULL,*ipointer=NULL,
66 *nx=NULL,*ny=NULL,*nz=NULL,symmetryflag=0,inputformat=0,mei[4],ido,ldz,
67 iparam[11],mxiter,info,ncv,jrow,lworkl,nev,ipntr[14],m,
68 *select=NULL,rvec=1,i,k,kflag,idesvar,node,inorm=0,irand=1,
69 iinc=1,noddiam=-1,ngraph,iobject,icoordinate,nrhs=1,ithermal=0,
70 ishape=0,inode,*irowd=NULL,*icold=NULL,*jqd=NULL,nzsc,*irowc=NULL,
71 *icolc=NULL,*jqc=NULL,nevold=0,imodes=0,niter=0;
72
73 double *xo=NULL,*yo=NULL,*zo=NULL,*x=NULL,*y=NULL,*z=NULL,*au=NULL,*ad=NULL,
74 *adb=NULL,*aub=NULL,*resid=NULL,*workd=NULL,*workl=NULL,tol,
75 *temp_array=NULL,*randval=NULL,fmin,fmax,pi,trace=0.0,dsum=0.,
76 *d=NULL,*zz=NULL,*acvector=NULL,*randfield=NULL,abserr,
77 actreliability,*stn=NULL,reliability,corrlen,*acscalar=NULL,
78 *add=NULL,*aud=NULL,sigma=0,*rhs=NULL,*vector=NULL,*adbd=NULL,
79 *aubd=NULL,*auc=NULL,delta,time=0.;
80
81 reliability=physcon[10];
82 corrlen=physcon[11];
83
84 /* Assigning values to eigenvalue parameters mei */
85
86 mei[1]=*ndesi; //number of requested Lanczos vectors
87 mei[2]=1000; //number of max iterations
88 mei[3]=0; //matrix storage
89
90 /* copying the frequency parameters */
91
92 pi=4.*atan(1.);
93 mxiter=1000;
94 tol=0.01;
95 fmin=2*pi*-1;
96 fmax=2*pi*-1;
97
98 /* prepare for near3d_se; routine "prefilter" can w/o any changes be taken
99 from the routine "filtermain.c" */
100
101 NNEW(xo,double,*ndesi);
102 NNEW(yo,double,*ndesi);
103 NNEW(zo,double,*ndesi);
104 NNEW(x,double,*ndesi);
105 NNEW(y,double,*ndesi);
106 NNEW(z,double,*ndesi);
107 NNEW(nx,ITG,*ndesi);
108 NNEW(ny,ITG,*ndesi);
109 NNEW(nz,ITG,*ndesi);
110
111 printf("\n Sorting the design variables for the computation of the covariance matrix\n\n");
112
113 for(m=0;m<*ndesi;m++){
114 xo[m]=co[3*(nodedesi[m]-1)];
115 x[m]=xo[m];
116 nx[m]=m+1;
117 yo[m]=co[3*(nodedesi[m]-1)+1];
118 y[m]=yo[m];
119 ny[m]=m+1;
120 zo[m]=co[3*(nodedesi[m]-1)+2];
121 z[m]=zo[m];
122 nz[m]=m+1;
123 }
124 kflag=2;
125 FORTRAN(dsort,(x,nx,ndesi,&kflag));
126 FORTRAN(dsort,(y,ny,ndesi,&kflag));
127 FORTRAN(dsort,(z,nz,ndesi,&kflag));
128
129 /* determining the structure of the covariance matrix */
130
131 printf(" Determining the structure of the covariance matrix\n\n");
132
133 nzss=20000000;
134 NNEW(mast1,ITG,nzss);
135 NNEW(irows,ITG,1);
136 NNEW(icols,ITG,*ndesi);
137 NNEW(jqs,ITG,*ndesi+1);
138 NNEW(ipointer,ITG,*ndesi);
139
140 mastructrand(icols,jqs,&mast1,&irows,ipointer,&nzss,ndesi,&corrlen,
141 xo,yo,zo,x,y,z,nx,ny,nz);
142
143 SFREE(mast1);SFREE(ipointer);
144 RENEW(irows,ITG,nzss);
145
146 SFREE(xo);SFREE(yo);SFREE(zo);SFREE(x);SFREE(y);SFREE(z);SFREE(nx);
147 SFREE(ny);SFREE(nz);
148
149 /* determining the entries of the autocovariance matrix */
150
151 printf(" Calculating the entries of the covariance matrix\n\n");
152
153 NNEW(ad,double,*ndesi);
154 NNEW(au,double,nzss);
155
156 FORTRAN(autocovmatrix,(co,ad,au,jqs,irows,ndesi,nodedesi,&corrlen,
157 randomval,irobustdesign));
158
159
160 /* ***************************************************************** */
161 /* start of conditional random field computation */
162 /* ***************************************************************** */
163
164 /* in case of a conditional random field, assemble the D-matrix and
165 the C-matrix */
166 /* have to be assembled: */
167 /* Cov^c = Cov - C * D^(-1) * C^(T) */
168 /* D-matrix: contains the covariance information of the boundary nodes */
169 /* C-matrix: contains the linking information between the design variables and
170 the boundary nodes subset) */
171
172 if(irobustdesign[2]==1){
173 nzsc=20000000;
174 NNEW(mast1,ITG,nzsc);
175 NNEW(irowc,ITG,1);
176 NNEW(icolc,ITG,*ndesi);
177 NNEW(jqc,ITG,*ndesi+1);
178 NNEW(ipointer,ITG,*ndesi);
179
180 mastructcmatrix(icolc,jqc,&mast1,&irowc,ipointer,&nzsc,ndesibou,
181 nodedesibou,nodedesiinvbou,jqs,irows,icols,ndesi,nodedesi);
182
183 SFREE(mast1);SFREE(ipointer);
184 RENEW(irowc,ITG,nzsc);
185
186 NNEW(auc,double,nzsc);
187
188 FORTRAN(cmatrix,(ad,au,jqs,irows,icols,ndesi,nodedesi,auc,jqc,irowc,
189 nodedesibou));
190
191 nzsd=20000000;
192 NNEW(mast1,ITG,nzsd);
193 NNEW(irowd,ITG,1);
194 NNEW(icold,ITG,*ndesi);
195 NNEW(jqd,ITG,*ndesi+1);
196 NNEW(ipointer,ITG,*ndesi);
197
198 mastructdmatrix(icold,jqd,&mast1,&irowd,ipointer,&nzsd,ndesibou,
199 nodedesibou,nodedesiinvbou,jqs,irows,icols,ndesi,nodedesi);
200
201 SFREE(mast1);SFREE(ipointer);
202 RENEW(irowd,ITG,nzsd);
203
204 printf(" Calculating the entries of the conditional covariance matrix\n\n");
205
206 NNEW(add,double,*ndesibou);
207 NNEW(aud,double,nzsd);
208
209 FORTRAN(dmatrix,(ad,au,jqs,irows,icols,ndesi,nodedesi,add,aud,
210 jqd,irowd,ndesibou,nodedesibou));
211
212 NNEW(adbd,double,*ndesibou);
213 NNEW(aubd,double,nzsd);
214
215 for(i=0;i<*ndesibou;i++){
216 adbd[i]=1.0;
217 }
218
219 /* LU decomposition of the left hand matrix */
220
221 if(*isolver==0){
222 #ifdef SPOOLES
223 spooles_factor(add,aud,adbd,aubd,&sigma,icold,irowd,ndesibou,&nzsd,
224 &symmetryflag,&inputformat,&nzsd);
225 #else
226 printf("*ERROR in randomfieldmain: the SPOOLES library is not linked\n\n");
227 FORTRAN(stop,());
228 #endif
229 }
230 else if(*isolver==4){
231 #ifdef SGI
232 token=1;
233 sgi_factor(add,aud,adbd,aubd,&sigma,icold,irowd,ndesibou,&nzsd,token);
234 #else
235 printf("*ERROR in randomfieldmain: the SGI library is not linked\n\n");
236 FORTRAN(stop,());
237 #endif
238 }
239 else if(*isolver==5){
240 #ifdef TAUCS
241 tau_factor(add,&aud,adbd,aubd,&sigma,icold,&irowd,ndesibou,&nzsd);
242 #else
243 printf("*ERROR in randomfieldmain: the TAUCS library is not linked\n\n");
244 FORTRAN(stop,());
245 #endif
246 }
247 else if(*isolver==7){
248 #ifdef PARDISO
249 pardiso_factor(add,aud,adbd,aubd,&sigma,icold,irowd,ndesibou,&nzsd,
250 &symmetryflag,&inputformat,jqd,&nzsd);
251 #else
252 printf("*ERROR in randomfieldmain: the PARDISO library is not linked\n\n");
253 FORTRAN(stop,());
254 #endif
255 }
256 else if(*isolver==8){
257 #ifdef PASTIX
258 pastix_factor_main(add,aud,adbd,aubd,&sigma,icold,irowd,ndesibou,&nzsd,
259 &symmetryflag,&inputformat,jqd,&nzsd);
260 #else
261 printf("*ERROR in randomfieldmain: the PASTIX library is not linked\n\n");
262 FORTRAN(stop,());
263 #endif
264 }
265
266 NNEW(rhs,double,*ndesibou);
267 NNEW(vector,double,*ndesi);
268
269 for(idesvar=1;idesvar<=*ndesi;idesvar++){
270
271 /* initialize rhs */
272
273 for(i=0;i<*ndesibou;i++){
274 rhs[i]=0.0;
275 }
276 for(i=0;i<*ndesi;i++){
277 vector[i]=0.0;
278 }
279
280 /* assembly of rhs */
281
282 FORTRAN(precondrandomfield,(auc,jqc,irowc,rhs,&idesvar));
283
284 /* solve the system */
285
286 if(*isolver==0){
287 #ifdef SPOOLES
288 spooles_solve(rhs,ndesibou);
289 #endif
290 }
291 else if(*isolver==4){
292 #ifdef SGI
293 sgi_solve(rhs,token);
294 #endif
295 }
296 else if(*isolver==5){
297 #ifdef TAUCS
298 tau_solve(rhs,ndesibou);
299 #endif
300 }
301 else if(*isolver==7){
302 #ifdef PARDISO
303 pardiso_solve(rhs,ndesibou,&symmetryflag,&inputformat,&nrhs);
304 #endif
305 }
306 else if(*isolver==8){
307 #ifdef PASTIX
308 pastix_solve(rhs,ndesibou,&symmetryflag,&nrhs);
309 #endif
310 }
311
312 /* carry out matrix multiplications */
313
314 FORTRAN(condrandomfield,(ad,au,jqs,irows,ndesi,rhs,vector,&idesvar,jqc,
315 auc,irowc));
316
317 }
318
319 /* clean the system */
320
321 if(*isolver==0){
322 #ifdef SPOOLES
323 spooles_cleanup();
324 #endif
325 }
326 else if(*isolver==4){
327 #ifdef SGI
328 sgi_cleanup(token);
329 #endif
330 }
331 else if(*isolver==5){
332 #ifdef TAUCS
333 tau_cleanup();
334 #endif
335 }
336 else if(*isolver==7){
337 #ifdef PARDISO
338 pardiso_cleanup(ndesibou,&symmetryflag,&inputformat);
339 #endif
340 }
341 else if(*isolver==8){
342 #ifdef PASTIX
343 #endif
344 }
345
346 SFREE(aubd);SFREE(adbd);SFREE(irowd);SFREE(jqd);SFREE(rhs);
347 SFREE(vector);SFREE(icold);SFREE(aud);SFREE(add);
348 SFREE(irowc);SFREE(jqc);SFREE(auc);SFREE(icolc);
349
350 }
351
352 /* ************************************************************ */
353 /* end of conditional random field computation */
354 /* ************************************************************ */
355
356
357 /* assign a starting value sigma for the calculation of the eigenvalues
358 sigma is set to trace(autocovmatrix) = maximum eigenvalue */
359
360 for(i=0;i<*ndesi;i++){
361 trace+=ad[i];
362 }
363
364 printf(" Trace of the covariance matrix: %e\n\n", trace);
365 /* Calculate the eigenmodes and eigenvalues of the (conditional)
366 autocovariance matrix */
367
368 NNEW(adb,double,*ndesi);
369 NNEW(aub,double,nzss);
370
371 for(i=0;i<*ndesi;i++){
372 adb[i]=1.0;
373 }
374
375 /* LU decomposition of the left hand matrix */
376
377 printf(" LU decomposition of the covariance matrix\n\n");
378
379 if(*isolver==0){
380 #ifdef SPOOLES
381 spooles_factor(ad,au,adb,aub,&trace,icols,irows,ndesi,&nzss,
382 &symmetryflag,&inputformat,&nzss);
383 #else
384 printf("*ERROR in randomfieldmain: the SPOOLES library is not linked\n\n");
385 FORTRAN(stop,());
386 #endif
387 }
388 else if(*isolver==4){
389 #ifdef SGI
390 token=1;
391 sgi_factor(ad,au,adb,aub,&trace,icols,irows,ndesi,&nzss,token);
392 #else
393 printf("*ERROR in randomfieldmain: the SGI library is not linked\n\n");
394 FORTRAN(stop,());
395 #endif
396 }
397 else if(*isolver==5){
398 #ifdef TAUCS
399 tau_factor(ad,&au,adb,aub,&trace,icols,&irows,ndesi,&nzss);
400 #else
401 printf("*ERROR in randomfieldmain: the TAUCS library is not linked\n\n");
402 FORTRAN(stop,());
403 #endif
404 }
405 else if(*isolver==6){
406 #ifdef MATRIXSTORAGE
407 matrixstorage(ad,&au,adb,aub,&trace,icols,&irows,ndesi,&nzss,
408 ntrans,inotr,trab,co,nk,nactdof,jobnamec,mi,ipkon,
409 lakon,kon,ne,mei,nboun,nmpc,cs,mcs,&ithermal,nmethod);
410 #else
411 printf("*ERROR in randomfieldmain: the MATRIXSTORAGE library is not linked\n\n");
412 FORTRAN(stop,());
413 #endif
414 }
415 else if(*isolver==7){
416 #ifdef PARDISO
417 pardiso_factor(ad,au,adb,aub,&trace,icols,irows,ndesi,&nzss,
418 &symmetryflag,&inputformat,jqs,&nzss);
419 #else
420 printf("*ERROR in randomfieldmain: the PARDISO library is not linked\n\n");
421 FORTRAN(stop,());
422 #endif
423 }
424 else if(*isolver==8){
425 #ifdef PASTIX
426 pastix_factor_main(ad,au,adb,aub,&trace,icols,irows,ndesi,&nzss,
427 &symmetryflag,&inputformat,jqs,&nzss);
428 #else
429 printf("*ERROR in randomfieldmain: the PASTIX library is not linked\n\n");
430 FORTRAN(stop,());
431 #endif
432 }
433
434 /* calculating the eigenvalues and eigenmodes */
435
436 /* initialization of relative error measure and number of eigenvalues*/
437
438 actreliability=0.;
439 nev=1;
440 ncv=4*nev;
441
442 NNEW(resid,double,*ndesi);
443 NNEW(zz,double,(long long)ncv**ndesi);
444 NNEW(workd,double,3**ndesi);
445 NNEW(workl,double,48);
446 NNEW(temp_array,double,*ndesi);
447 NNEW(select,ITG,4);
448 NNEW(d,double,nev);
449
450 while(actreliability<=reliability){
451
452 niter+=1;
453
454 if(nev>=*ndesi){
455 printf("*WARNING in randomfieldmain: number of requested\n");
456 printf(" eigenvectors must be smaller than the \n");
457 printf(" number of design variables; \n");
458 printf(" the number of eigenvectors is set to \n");
459 printf(" the number of design variables minus one. \n");
460 nev=*ndesi-1;
461 break;
462 }
463
464 printf(" Calculating the eigenvalues and the eigenmodes\n\n");
465
466 ncv=4*nev;
467 ido=0;
468 ldz=*ndesi;
469 iparam[0]=1;
470 iparam[2]=mxiter;
471 iparam[3]=1;
472 iparam[6]=3;
473 info=0;
474
475 RENEW(zz,double,(long long)ncv**ndesi);
476
477 lworkl=ncv*(8+ncv);
478 RENEW(workl,double,lworkl);
479 FORTRAN(dsaupd,(&ido,bmat,ndesi,which,&nev,&tol,resid,&ncv,zz,&ldz,
480 iparam,ipntr,workd,workl,&lworkl,&info));
481
482 while((ido==-1)||(ido==1)||(ido==2)){
483 if(ido==-1){
484 FORTRAN(op,(ndesi,&workd[ipntr[0]-1],temp_array,adb,aub,jqs,irows));
485 }
486 if((ido==-1)||(ido==1)){
487
488 /* solve the linear equation system */
489
490 if(ido==-1){
491 if(*isolver==0){
492 #ifdef SPOOLES
493 spooles_solve(temp_array,ndesi);
494 #endif
495 }
496 else if(*isolver==4){
497 #ifdef SGI
498 sgi_solve(temp_array,token);
499 #endif
500 }
501 else if(*isolver==5){
502 #ifdef TAUCS
503 tau_solve(temp_array,ndesi);
504 #endif
505 }
506 else if(*isolver==7){
507 #ifdef PARDISO
508 pardiso_solve(temp_array,ndesi,&symmetryflag,&inputformat,&nrhs);
509 #endif
510 }
511 else if(*isolver==8){
512 #ifdef PASTIX
513 pastix_solve(temp_array,ndesi,&symmetryflag,&nrhs);
514 #endif
515 }
516
517 for(jrow=0;jrow<*ndesi;jrow++){
518 workd[ipntr[1]-1+jrow]=temp_array[jrow];
519 }
520 }
521 else if(ido==1){
522 if(*isolver==0){
523 #ifdef SPOOLES
524 spooles_solve(&workd[ipntr[2]-1],ndesi);
525 #endif
526 }
527 else if(*isolver==4){
528 #ifdef SGI
529 sgi_solve(&workd[ipntr[2]-1],token);
530 #endif
531 }
532 else if(*isolver==5){
533 #ifdef TAUCS
534 tau_solve(&workd[ipntr[2]-1],ndesi);
535 #endif
536 }
537 else if(*isolver==7){
538 #ifdef PARDISO
539 pardiso_solve(&workd[ipntr[2]-1],ndesi,&symmetryflag,&inputformat,&nrhs);
540 #endif
541 }
542 else if(*isolver==8){
543 #ifdef PASTIX
544 pastix_solve(&workd[ipntr[2]-1],ndesi,&symmetryflag,&nrhs);
545 #endif
546 }
547
548 for(jrow=0;jrow<*ndesi;jrow++){
549 workd[ipntr[1]-1+jrow]=workd[ipntr[2]-1+jrow];
550 }
551 }
552
553 }
554
555 if(ido==2){
556 FORTRAN(op,(ndesi,&workd[ipntr[0]-1],&workd[ipntr[1]-1],
557 adb,aub,jqs,irows));
558 }
559
560 FORTRAN(dsaupd,(&ido,bmat,ndesi,which,&nev,&tol,resid,&ncv,zz,&ldz,
561 iparam,ipntr,workd,workl,&lworkl,&info));
562 }
563
564 RENEW(select,ITG,ncv);
565 RENEW(d,double,nev);
566
567 FORTRAN(dseupd,(&rvec,howmny,select,d,zz,&ldz,&trace,bmat,ndesi,which,
568 &nev,&tol,resid,&ncv,zz,&ldz,iparam,ipntr,workd,workl,
569 &lworkl,&info));
570
571 /* calculating the error measure and estimation of number of
572 eigenvalues for next iteration */
573
574 for(i=nev-nevold-1;i>=0;i--){
575 imodes+=1;
576 dsum+=d[i];
577 actreliability=dsum/trace;
578 FORTRAN(writerandomfield,(&d[i],&actreliability,&imodes));
579 }
580
581 printf("\n Reliability iteration: %" ITGFORMAT "\n", niter);
582 printf(" Number of eigenvalues calculated: %" ITGFORMAT "\n", nev);
583 printf(" Corresponding reliability: %e\n\n", actreliability);
584
585 nevold=nev;
586 delta=reliability-actreliability;
587 nev=(ITG)(nev+delta/(d[0]/trace)+1);
588
589 } /*while loop end*/
590
591 /* storing the mean in the frd file */
592
593 NNEW(acvector,double,3**nk);
594 NNEW(acscalar,double,*nk);
595
596 ++*kode;
597 for(idesvar=0;idesvar<*ndesi;idesvar++){
598 node=nodedesi[idesvar]-1;
599 if(irobustdesign[1]==1){
600 acscalar[node]=randomval[0];
601 for(i=0;i<3;i++){
602 acvector[3*node+i]=randomval[0]*extnor[3*node+i];
603 }
604 }else{
605 acscalar[node]=randomval[2*node];
606 for(i=0;i<3;i++){
607 acvector[3*node+i]=randomval[2*node]*extnor[3*node+i];
608 }
609 }
610 }
611
612 irand=3;
613 frd_sen(co,nk,stn,inum,nmethod,kode,filab,&time,nstate_,istep,
614 &iinc,&k,&noddiam,description,mi,&ngraph,ne,cs,set,nset,
615 istartset,iendset,ialset,jobnamec,output,
616 acscalar,&iobject,objectset,ntrans,inotr,trab,&idesvar,orname,
617 &icoordinate,&inorm,&irand,&ishape);
618 irand=4;
619 frd_sen(co,nk,stn,inum,nmethod,kode,filab,&time,nstate_,istep,
620 &iinc,&k,&noddiam,description,mi,&ngraph,ne,cs,set,nset,
621 istartset,iendset,ialset,jobnamec,output,
622 acvector,&iobject,objectset,ntrans,inotr,trab,&idesvar,orname,
623 &icoordinate,&inorm,&irand,&ishape);
624 irand=1;
625
626 /* storing the eigenvectors of the autocovariance matrix in the frd file
627 acvector = autocovariance vector */
628
629 for(k=nev;k>=1;k--){
630
631 /* generating and writing the autocorrelation vector */
632
633 ++*kode;
634 for(idesvar=0;idesvar<*ndesi;idesvar++){
635 node=nodedesi[idesvar]-1;
636 acscalar[node]=zz[(k-1)**ndesi+idesvar];
637 for(i=0;i<3;i++){
638 acvector[3*node+i]=zz[(k-1)**ndesi+idesvar]*extnor[3*node+i];
639 }
640 }
641
642 frd_sen(co,nk,stn,inum,nmethod,kode,filab,&d[k-1],nstate_,
643 istep,
644 &iinc,&k,&noddiam,description,mi,&ngraph,ne,cs,set,nset,
645 istartset,iendset,ialset,jobnamec,output,
646 acscalar,&iobject,objectset,ntrans,inotr,trab,&idesvar,orname,
647 &icoordinate,&inorm,&irand,&ishape);
648 irand=2;
649 frd_sen(co,nk,stn,inum,nmethod,kode,filab,&d[k-1],nstate_,
650 istep,
651 &iinc,&k,&noddiam,description,mi,&ngraph,ne,cs,set,nset,
652 istartset,iendset,ialset,jobnamec,output,
653 acvector,&iobject,objectset,ntrans,inotr,trab,&idesvar,orname,
654 &icoordinate,&inorm,&irand,&ishape);
655 irand=1;
656
657 }
658
659 /* generating and writing the local reliability */
660
661 ++*kode;
662 for(idesvar=0;idesvar<*ndesi;idesvar++){
663 node=nodedesi[idesvar]-1;
664 acscalar[node]=0.0;
665 }
666 for(k=nev;k>=1;k--){
667 for(idesvar=0;idesvar<*ndesi;idesvar++){
668 node=nodedesi[idesvar]-1;
669 acscalar[node]+=zz[(k-1)**ndesi+idesvar]*zz[(k-1)**ndesi+idesvar]*
670 d[k-1];
671 }
672 }
673
674 for(idesvar=0;idesvar<*ndesi;idesvar++){
675 node=nodedesi[idesvar]-1;
676 if(irobustdesign[1]==1){
677 acscalar[node]=acscalar[node]/(randomval[1]*randomval[1]);
678 }else{
679 acscalar[node]=acscalar[node]/(randomval[2*(node)+1]*
680 randomval[2*(node)+1]);
681 }
682 }
683
684 irand=5;
685 frd_sen(co,nk,stn,inum,nmethod,kode,filab,&time,nstate_,
686 istep,
687 &iinc,&k,&noddiam,description,mi,&ngraph,ne,cs,set,nset,
688 istartset,iendset,ialset,jobnamec,output,
689 acscalar,&iobject,objectset,ntrans,inotr,trab,&idesvar,orname,
690 &icoordinate,&inorm,&irand,&ishape);
691 /*
692 -----------
693 free memory
694 -----------
695 */
696
697 SFREE(temp_array);
698 if(*isolver==0){
699 #ifdef SPOOLES
700 spooles_cleanup();
701 #endif
702 }
703 else if(*isolver==4){
704 #ifdef SGI
705 sgi_cleanup(token);
706 #endif
707 }
708 else if(*isolver==5){
709 #ifdef TAUCS
710 tau_cleanup();
711 #endif
712 }
713 else if(*isolver==7){
714 #ifdef PARDISO
715 pardiso_cleanup(ndesi,&symmetryflag,&inputformat);
716 #endif
717 }
718 else if(*isolver==8){
719 #ifdef PASTIX
720 #endif
721 }
722
723 if(info!=0){
724 printf("*ERROR in randomfieldmain: info=%" ITGFORMAT "\n",info);
725 printf(" # of converged eigenvalues=%" ITGFORMAT "\n\n",iparam[4]);
726 }
727
728 SFREE(select);SFREE(workd);SFREE(workl);SFREE(resid);
729 SFREE(irows);SFREE(icols);SFREE(jqs);SFREE(au);SFREE(ad);
730 SFREE(adb);SFREE(aub);SFREE(d);SFREE(zz);SFREE(randval);
731 SFREE(randfield);SFREE(acvector);SFREE(acscalar);
732
733 return;
734
735 }
736
737 #endif
738