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 <stdio.h>
19 #include <math.h>
20 #include <stdlib.h>
21 #include "CalculiX.h"
22
preiter(double * ad,double ** aup,double * b,ITG ** icolp,ITG ** irowp,ITG * neq,ITG * nzs,ITG * isolver,ITG * iperturb)23 void preiter(double *ad, double **aup, double *b, ITG **icolp, ITG **irowp,
24 ITG *neq, ITG *nzs, ITG *isolver, ITG *iperturb){
25
26 ITG precFlg,niter=5000000,ndim,i,j,k,ier,*icol=NULL,*irow=NULL,
27 *irow_save=NULL,*icol_save=NULL;
28 double eps=1.e-4,*u=NULL,*au=NULL;
29
30 if(*neq==0) return;
31
32 /* icol(i) = # subdiagonal nonzeros in column i (i=1,neq)
33 irow(i) = row number of entry i in au (i=1,nzs)
34 ad(i) = diagonal term in column i of the matrix
35 au(i) = subdiagonal nonzero term i; the terms are entered
36 column per column */
37
38 au=*aup;
39 irow=*irowp;
40 icol=*icolp;
41
42 if(*iperturb>1){
43 NNEW(irow_save,ITG,*nzs);
44 NNEW(icol_save,ITG,*neq);
45 for(i=0;i<*nzs;++i){
46 irow_save[i]=irow[i];
47 }
48 for(i=0;i<*neq;++i){
49 icol_save[i]=icol[i];
50 }
51 }
52
53 if(*isolver==2) {precFlg=0;}
54 else {precFlg=3;}
55
56 ndim=*neq+*nzs;
57
58 RENEW(au,double,ndim);
59 RENEW(irow,ITG,ndim);
60 RENEW(icol,ITG,ndim);
61
62 k=*nzs;
63 for(i=*neq-1;i>=0;--i){
64 for(j=0;j<icol[i];++j){
65 icol[--k]=i+1;
66 }
67 }
68
69 k=*nzs;
70 j=0;
71 for(i=0;i<*neq;++i){
72 au[k]=ad[i];
73 irow[k]=++j;
74 icol[k]=j;
75 ++k;
76 }
77
78 /* rearranging the storage of the left hand side */
79
80 FORTRAN(rearrange,(au,irow,icol,&ndim,neq));
81
82 RENEW(irow,ITG,*neq);
83
84 NNEW(u,double,*neq);
85
86 ier=cgsolver(au,u,b,*neq,ndim,icol,irow,&eps,&niter,precFlg);
87
88 printf("error condition (0=good, 1=bad) = %" ITGFORMAT "\n",ier);
89 printf("# of iterations = %" ITGFORMAT "\n",niter);
90
91 for(i=0;i<*neq;++i){
92 b[i]=u[i];
93 }
94
95 SFREE(u);
96
97 if(*iperturb>1){
98 RENEW(irow,ITG,*nzs);
99 RENEW(icol,ITG,*neq);
100 for(i=0;i<*nzs;++i){
101 irow[i]=irow_save[i];
102 }
103 for(i=0;i<*neq;++i){
104 icol[i]=icol_save[i];
105 }
106 SFREE(irow_save);SFREE(icol_save);
107 }
108
109 *aup=au;
110 *irowp=irow;
111 *icolp=icol;
112
113 return;
114 }
115