CalculiX  2.8
A Free Software Three-Dimensional Structural Finite Element Program
 All Classes Files Functions Variables Macros
pcgsolver.c File Reference
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "CalculiX.h"
Include dependency graph for pcgsolver.c:

Go to the source code of this file.

Macros

#define GOOD   0
 
#define BAD   1
 
#define FALSE   0
 
#define TRUE   1
 

Functions

ITG cgsolver (double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, ITG *iz, double *eps, ITG *niter, ITG precFlg)
 
void PCG (double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, ITG *iz, double *eps, ITG *niter, ITG precFlg, double *rho, double *r, double *g, double *C, double *z)
 
void CG (double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, ITG *iz, double *eps, ITG *niter, double *r, double *p, double *z)
 
void Scaling (double *A, double *b, ITG neq, ITG *ia, ITG *iz, double *d)
 
void MatVecProduct (double *A, double *p, ITG neq, ITG *ia, ITG *iz, double *z)
 
void PreConditioning (double *A, ITG neq, ITG len, ITG *ia, ITG *iz, double alpha, ITG precFlg, double *C, ITG *ier)
 
void Mrhor (double *C, ITG neq, ITG *ia, ITG *iz, double *r, double *rho)
 
void InnerProduct (double *a, double *b, ITG n, double *Sum)
 

Macro Definition Documentation

#define BAD   1

Definition at line 26 of file pcgsolver.c.

#define FALSE   0

Definition at line 27 of file pcgsolver.c.

#define GOOD   0

Definition at line 25 of file pcgsolver.c.

#define TRUE   1

Definition at line 28 of file pcgsolver.c.

Function Documentation

void CG ( double *  A,
double *  x,
double *  b,
ITG  neq,
ITG  len,
ITG ia,
ITG iz,
double *  eps,
ITG niter,
double *  r,
double *  p,
double *  z 
)

Definition at line 504 of file pcgsolver.c.

506 {
507  ITG i=0, k=0, ncg=0,iam;
508  double ekm1=0.0,c1=0.005,qam,ram=0.,err;
509  double rr=0.0, pz=0.0, qk=0.0, rro=0.0;
510 
511 
512  /* solving the system of equations using the iterative solver */
513 
514  printf("Solving the system of equations using the iterative solver\n\n");
515 
516 /*..initialize result, search and residual vectors................................. */
517  qam=0.;iam=0;
518  for (i=0; i<neq; i++)
519  {
520  x[i] = 0.0; /*..start vector x0 = 0.................... */
521  r[i] = b[i]; /*..residual vector r0 = Ax+b = b.......... */
522  p[i] = -r[i]; /*..initial search vector.................. */
523  err=fabs(r[i]);
524  if(err>1.e-20){qam+=err;iam++;}
525  }
526  if(iam>0) qam=qam/neq;
527  else {*niter=0;return;}
528  /*else qam=0.01;*/
529 /*..main iteration loop............................................................ */
530  for (k=1; k<=(*niter); k++, ncg++)
531  {
532 /*......inner product rT r......................................................... */
533  InnerProduct(r,r,neq,&rr);
534  printf("iteration= %" ITGFORMAT ", error= %e, limit=%e\n",ncg,ram,c1*qam);
535 /*......If working with Penalty-terms for boundary conditions you can get nume-.... */
536 /*......rical troubles because RR becomes a very large value. With at least two.... */
537 /*......iterations the result may be better !!!.................................... */
538 /*......convergency check.......................................................... */
539  if (k!=1 && (ram<=c1*qam)) break;
540 /*......new search vector.......................................................... */
541  if (k!=1)
542  {
543  ekm1 = rr/rro;
544  for (i=0; i<neq; i++) p[i] = ekm1*p[i]-r[i];
545  }
546 /*......matrix vector product A p = z.............................................. */
547  MatVecProduct(A,p,neq,ia,iz,z);
548 /*......inner product pT z......................................................... */
549  InnerProduct(p,z,neq,&pz);
550 /*......step size.................................................................. */
551  qk = rr/pz;
552 /*......approximated solution vector, residual vector.............................. */
553  ram=0.;
554  for (i=0; i<neq; i++)
555  {
556  x[i] = x[i] + qk*p[i];
557  r[i] = r[i] + qk*z[i];
558  err=fabs(r[i]);
559  if(err>ram) ram=err;
560  }
561 /*......store previous residual.................................................... */
562  rro = rr;
563 
564  }
565  if(k==*niter){
566  printf("*ERROR in PCG: no convergence;");
567  FORTRAN(stop,());
568  }
569  *eps = rr; /*..return residual............................ */
570  *niter = ncg; /*..return number of iterations................ */
571 /*..That's it...................................................................... */
572  return;
573 }
#define ITGFORMAT
Definition: CalculiX.h:52
void FORTRAN(addimdnodecload,(ITG *nodeforc, ITG *i, ITG *imdnode, ITG *nmdnode, double *xforc, ITG *ikmpc, ITG *ilmpc, ITG *ipompc, ITG *nodempc, ITG *nmpc, ITG *imddof, ITG *nmddof, ITG *nactdof, ITG *mi, ITG *imdmpc, ITG *nmdmpc, ITG *imdboun, ITG *nmdboun, ITG *ikboun, ITG *nboun, ITG *ilboun, ITG *ithermal))
void MatVecProduct(double *A, double *p, ITG neq, ITG *ia, ITG *iz, double *z)
Definition: pcgsolver.c:323
subroutine stop()
Definition: stop.f:19
void InnerProduct(double *a, double *b, ITG n, double *Sum)
Definition: pcgsolver.c:477
#define ITG
Definition: CalculiX.h:51
ITG cgsolver ( double *  A,
double *  x,
double *  b,
ITG  neq,
ITG  len,
ITG ia,
ITG iz,
double *  eps,
ITG niter,
ITG  precFlg 
)

Definition at line 70 of file pcgsolver.c.

73 {
74  ITG i=0;
75  double *Factor=NULL,*r=NULL,*p=NULL,*z=NULL,*C=NULL,*g=NULL,*rho=NULL;
76 
77  /* reduce row and column indices by 1 (FORTRAN->C) */
78 
79  for (i=0; i<neq; i++) --iz[i];
80  for (i=0; i<len; i++) --ia[i];
81 
82  /* Scaling the equation system A x + b = 0 */
83 
84  NNEW(Factor,double,neq);
85  Scaling(A,b,neq,ia,iz,Factor);
86 
87  /* SOLVER/PRECONDITIONING TYPE */
88 
89  /* Conjugate gradient solver without preconditioning */
90 
91  if (!precFlg)
92  {
93  NNEW(r,double,neq);
94  NNEW(p,double,neq);
95  NNEW(z,double,neq);
96  CG(A,x,b,neq,len,ia,iz,eps,niter,r,p,z);
97  SFREE(r);SFREE(p);SFREE(z);
98  }
99 
100  /* Conjugate gradient solver with incomplete Cholesky preconditioning on
101  full matrix */
102 
103  else if (precFlg==3)
104  {
105  NNEW(rho,double,neq);
106  NNEW(r,double,neq);
107  NNEW(g,double,neq);
108  NNEW(C,double,len);
109  NNEW(z,double,neq);
110  PCG(A,x,b,neq,len,ia,iz,eps,niter,precFlg,rho,r,g,C,z);
111  SFREE(rho);SFREE(r);SFREE(g);SFREE(C);SFREE(z);
112  }
113 
114  /* Backscaling of the solution vector */
115 
116  for (i=0; i<neq; i++) x[i] *= Factor[i];
117 
118  /* That's it */
119 
120  SFREE(Factor);
121  return GOOD;
122 }
#define GOOD
Definition: pcgsolver.c:25
void CG(double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, ITG *iz, double *eps, ITG *niter, double *r, double *p, double *z)
Definition: pcgsolver.c:504
void Scaling(double *A, double *b, ITG neq, ITG *ia, ITG *iz, double *d)
Definition: pcgsolver.c:281
#define SFREE(a)
Definition: CalculiX.h:41
#define ITG
Definition: CalculiX.h:51
#define NNEW(a, b, c)
Definition: CalculiX.h:39
void PCG(double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, ITG *iz, double *eps, ITG *niter, ITG precFlg, double *rho, double *r, double *g, double *C, double *z)
Definition: pcgsolver.c:147
void InnerProduct ( double *  a,
double *  b,
ITG  n,
double *  Sum 
)

Definition at line 477 of file pcgsolver.c.

478 {
479  ITG i=0;
480 /*..local vectors.................................................................. */
481  *Sum=0.;
482  for (i=0; i<n; i++) *Sum += a[i]*b[i];
483  return;
484 }
#define ITG
Definition: CalculiX.h:51
void MatVecProduct ( double *  A,
double *  p,
ITG  neq,
ITG ia,
ITG iz,
double *  z 
)

Definition at line 323 of file pcgsolver.c.

325 {
326  ITG i=0, j=0, jlo=0, jup=0, k=0;
327 
328  /* matrix vector product */
329 
330  z[0] = A[iz[0]]*p[0];
331  for (i=1; i<neq; i++)
332  {
333  z[i] = A[iz[i]]*p[i];
334 
335  /* first non-zero element in current row */
336 
337  jlo = iz[i-1]+1;
338 
339  /* last non-zero off-diagonal element */
340 
341  jup = iz[i]-1;
342  for (j=jlo; j<=jup; j++)
343  {
344  k = ia[j];
345  z[i] += A[j]*p[k];
346  z[k] += A[j]*p[i];
347  }
348  }
349  return;
350 }
#define ITG
Definition: CalculiX.h:51
void Mrhor ( double *  C,
ITG  neq,
ITG ia,
ITG iz,
double *  r,
double *  rho 
)

Definition at line 445 of file pcgsolver.c.

446 {
447  ITG i=0, j=0, jlo=0, jup=0;
448  double s=0.0;
449 /*..solve equation sytem by forward/backward substitution.......................... */
450  rho[0] = r[0];
451  for (i=1; i<neq; i++)
452  {
453  s = 0.0;
454  jlo = iz[i-1]+1; /*..first non-zero element in current row...... */
455  jup = iz[i]; /*..diagonal element in current row............ */
456  for (j=jlo; j<jup; j++) /*..all non-zero off-diagonal element.......... */
457  s += C[j]*rho[ia[j]];
458  rho[i] = (r[i]-s)/C[jup];
459  }
460  for (i=neq-1; i>0; i--)
461  {
462  rho[i] /= C[iz[i]];
463  jlo = iz[i-1]+1; /*..first non-zero element in current row...... */
464  jup = iz[i]-1; /*..diagonal element in current row............ */
465  for (j=jlo; j<=jup; j++) /*..all non-zero off-diagonal element.......... */
466  rho[ia[j]] -= C[j]*rho[i];
467  }
468  return;
469 }
#define ITG
Definition: CalculiX.h:51
void PCG ( double *  A,
double *  x,
double *  b,
ITG  neq,
ITG  len,
ITG ia,
ITG iz,
double *  eps,
ITG niter,
ITG  precFlg,
double *  rho,
double *  r,
double *  g,
double *  C,
double *  z 
)

Definition at line 147 of file pcgsolver.c.

150 {
151  ITG i=0, k=0, ncg=0,iam,ier=0;
152  double alpha=0.0, ekm1=0.0, rrho=0.0;
153  double rrho1=0.0, gz=0.0, qk=0.0;
154  double c1=0.005,qam,err,ram=0;
155 
156 
157  /* initialize result and residual vectors */
158 
159  qam=0.;iam=0;
160  for (i=0; i<neq; i++)
161  {
162  x[i] = 0.0;
163  r[i] = b[i];
164  err=fabs(r[i]);
165  if(err>1.e-20){qam+=err;iam++;}
166  }
167  if(iam>0) qam=qam/iam;
168  else {*niter=0;return;}
169 
170  /* preconditioning */
171 
172  printf("Cholesky preconditioning\n\n");
173 
174  printf("alpha=%f\n",alpha);
175  PreConditioning(A,neq,len,ia,iz,alpha,precFlg,C,&ier);
176  while (ier==0)
177  {
178  if (alpha<=0.0) alpha=0.005;
179  alpha += alpha;
180  printf("alpha=%f\n",alpha);
181  PreConditioning(A,neq,len,ia,iz,alpha,precFlg,C,&ier);
182  }
183 
184  /* solving the system of equations using the iterative solver */
185 
186  printf("Solving the system of equations using the iterative solver\n\n");
187 
188  /* main iteration loop */
189 
190  for (k=1; k<=*niter; k++, ncg++)
191  {
192 
193  /* solve M rho = r, M=C CT */
194 
195  Mrhor(C,neq,ia,iz,r,rho);
196 
197  /* inner product (r,rho) */
198 
199  InnerProduct(r,rho,neq,&rrho);
200 
201  /* If working with Penalty-terms for boundary conditions you can get
202  numerical troubles because RRHO becomes a very large value.
203  With at least two iterations the result may be better !!! */
204 
205  /* convergency check */
206 
207  printf("iteration= %" ITGFORMAT ", error= %e, limit=%e\n",ncg,ram,c1*qam);
208  if (k!=1 && (ram<=c1*qam)) break;
209  if (k!=1)
210  {
211  ekm1 = rrho/rrho1;
212  for (i=0; i<neq; i++) g[i] = ekm1*g[i]-rho[i];
213  }
214  else
215  {
216 
217  /* g1 = -rho0 */
218 
219  for (i=0; i<neq; i++) g[i] = -rho[i];
220  }
221 
222  /* solve equation system z = A g_k */
223 
224  MatVecProduct(A,g,neq,ia,iz,z);
225 
226  /* inner product (g,z) */
227 
228  InnerProduct(g,z,neq,&gz);
229  qk = rrho/gz;
230  ram=0.;
231  for (i=0; i<neq; i++)
232  {
233  x[i] += qk*g[i];
234  r[i] += qk*z[i];
235  err=fabs(r[i]);
236  if(err>ram) ram=err;
237  }
238  rrho1 = rrho;
239  }
240  if(k==*niter){
241  printf("*ERROR in PCG: no convergence;");
242  FORTRAN(stop,());
243  }
244  *eps = rrho;
245  *niter = ncg;
246 
247  return;
248 }
#define ITGFORMAT
Definition: CalculiX.h:52
void FORTRAN(addimdnodecload,(ITG *nodeforc, ITG *i, ITG *imdnode, ITG *nmdnode, double *xforc, ITG *ikmpc, ITG *ilmpc, ITG *ipompc, ITG *nodempc, ITG *nmpc, ITG *imddof, ITG *nmddof, ITG *nactdof, ITG *mi, ITG *imdmpc, ITG *nmdmpc, ITG *imdboun, ITG *nmdboun, ITG *ikboun, ITG *nboun, ITG *ilboun, ITG *ithermal))
void MatVecProduct(double *A, double *p, ITG neq, ITG *ia, ITG *iz, double *z)
Definition: pcgsolver.c:323
subroutine stop()
Definition: stop.f:19
void InnerProduct(double *a, double *b, ITG n, double *Sum)
Definition: pcgsolver.c:477
void PreConditioning(double *A, ITG neq, ITG len, ITG *ia, ITG *iz, double alpha, ITG precFlg, double *C, ITG *ier)
Definition: pcgsolver.c:373
void Mrhor(double *C, ITG neq, ITG *ia, ITG *iz, double *r, double *rho)
Definition: pcgsolver.c:445
#define ITG
Definition: CalculiX.h:51
void PreConditioning ( double *  A,
ITG  neq,
ITG  len,
ITG ia,
ITG iz,
double  alpha,
ITG  precFlg,
double *  C,
ITG ier 
)

Definition at line 373 of file pcgsolver.c.

376 {
377  ITG i=0, j=0, jlo=0, jup=0, k=0, klo=0, kup=0, l=0, llo=0, lup=0;
378  ITG id=0, nILU=0, m=0;
379  double factor;
380  factor = 1.0/(1.0+alpha);
381 /*..division of the off-diagonal elements of A by factor (1.0+alpha)............... */
382  C[0] = A[0];
383  for (i=1; i<neq; i++)
384  {
385  jlo = iz[i-1]+1; /*..first non-zero element in current row...... */
386  jup = iz[i]; /*..diagonal element........................... */
387  C[jup] = A[jup]; /*..copy of diagonal element................... */
388  for (j=jlo; j<jup; j++) /*..all non-zero off-diagonal elements......... */
389  C[j] = A[j]*factor;
390  }
391  nILU = neq; /*..ILU on full matrix................. */
392 /*..partial Cholesky decomposition of scaled matrix A.............................. */
393  for (i=1; i<nILU; i++)
394  {
395  jlo = iz[i-1]+1; /*..first non-zero element in current row...... */
396  jup = iz[i]; /*..diagonal element........................... */
397  for (j=jlo; j<jup; j++) /*..all non-zero off-diagonal elements......... */
398  {
399  C[j] /= C[iz[ia[j]]];
400  klo = j+1; /*..next non-zero element in current row....... */
401  kup = jup; /*..diagonal element in current row............ */
402  for (k=klo; k<=kup; k++)
403  {
404  m = ia[k];
405  llo = iz[m-1]+1;
406  lup = iz[m];
407  for (l=llo; l<=lup; l++)
408  {
409  if (ia[l]>ia[j]) break;
410  if (ia[l]<ia[j]) continue;
411  C[k] -= C[j]*C[l];
412  break;
413  }
414  }
415  }
416  id = iz[i];
417  if (C[id]<1.0e-6){
418  return;
419  }
420  C[id] = sqrt(C[id]);
421  }
422  *ier=1;
423  return;
424 }
#define ITG
Definition: CalculiX.h:51
void Scaling ( double *  A,
double *  b,
ITG  neq,
ITG ia,
ITG iz,
double *  d 
)

Definition at line 281 of file pcgsolver.c.

282 {
283  ITG i=0, j=0, jlo=0, jup=0;
284 
285  /* extract diagonal vector from matrix A */
286 
287  for (i=0; i<neq; i++) d[i] = 1.0/sqrt(A[iz[i]]);
288 
289  /* scale right hand side (Ax=b -> Ax+b=0: negative sign) */
290 
291  for (i=0; i<neq; i++) b[i] *= -d[i];
292 
293  /* scale matrix A */
294 
295  A[iz[0]] *= d[0]*d[0];
296  for (i=1; i<neq; i++)
297  {
298  jlo = iz[i-1]+1;
299  jup = iz[i];
300  for (j=jlo; j<=jup; j++) A[j] *= d[i]*d[ia[j]];
301  }
302 
303  return;
304 }
#define ITG
Definition: CalculiX.h:51
Hosted by OpenAircraft.com, (Michigan UAV, LLC)