root/src/preiter.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
```   1 /*     CalculiX - A 3-dimensional finite element program                   */
2 /*              Copyright (C) 1998-2015 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
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){
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 }
```

/* [<][>][^][v][top][bottom][index][help] */