root/src/pcgsolver.c

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

DEFINITIONS

This source file includes following definitions.
  1. cgsolver
  2. PCG
  3. Scaling
  4. MatVecProduct
  5. PreConditioning
  6. Mrhor
  7. InnerProduct
  8. CG

   1 /*
   2 
   3   pcgsolver.c   --      The preconditioned conjugate gradient solver(s)
   4   
   5                                 Implemented solver/preconditioner:
   6                                     preconditioned cg-solver with
   7                                        * diagonal scaling only
   8                                        * incomplete Cholesky on full matrix
   9                                 
  10                                 Most functions are based on:
  11                                 H.R. Schwarz: FORTRAN-Programme zur Methode 
  12                                              der finiten Elemente,
  13                                               Teubner, 1981 
  14  
  15                                The present version is based on the c-code by
  16                                Martin Ruecker and Ernst Rank of the Technical
  17                                University of Munich 
  18  */
  19 
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23 #include "CalculiX.h"
  24 
  25 #define GOOD 0
  26 #define BAD 1
  27 #define FALSE 0
  28 #define TRUE 1
  29 
  30 /* Prototyping  */
  31 
  32 ITG cgsolver (double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, 
  33                ITG *iz,double *eps, ITG *niter, ITG precFlg);
  34 void PCG (double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, 
  35           ITG *iz,double *eps, ITG *niter, ITG precFlg,
  36           double *rho, double *r, double *g, double *C, double *z);
  37 void CG (double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, 
  38          ITG *iz,double *eps, ITG *niter,double *r, double *p, double *z);
  39 void Scaling (double *A, double *b, ITG neq, ITG *ia, ITG *iz, double *d);
  40 void MatVecProduct (double *A, double *p, ITG neq, ITG *ia, ITG *iz, 
  41                     double *z);
  42 void PreConditioning (double *A, ITG neq, ITG len, ITG *ia, ITG *iz, 
  43                       double alpha, ITG precFlg,double *C, ITG *ier);
  44 void Mrhor (double *C, ITG neq, ITG *ia, ITG *iz, double *r, double *rho);
  45 void InnerProduct (double *a, double *b, ITG n, double *Sum);
  46         
  47 /* **********************************************************************
  48 
  49 The (preconditioned) conjugate gradient solver
  50 
  51   parameter:                                                                
  52            A       compact row oriented storage of lower left of matrix A  
  53            neq     order of A, number of equations                         
  54            len     number of non zero entries in Matrix A                  
  55            ia      column indices of corresponding elements in Matrix A    
  56            iz      row indices (diagonal elements indices)               
  57            x       solution vector                                       
  58            b       right hand side                                       
  59            eps     required accuracy -> residual                         
  60            niter   maximum number of iterations -> number of iterations  
  61            precFlg preconditioning flag                                  
  62 
  63 The compact row oriented storage of sparse quadratic matrices is decsribed in
  64 H.R. Schwarz: FORTRAN-Programme zur Methode der finiten Elemente, pp.66-67, 
  65 Teubner, 1981
  66  
  67 ********************************************************************** 
  68 */
  69 
  70 ITG cgsolver (double *A, double *x, double *b, ITG neq, ITG len, 
  71                ITG *ia, ITG *iz, 
  72                                 double *eps, ITG *niter, ITG precFlg)
  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 }
 123 
 124 /* **********************************************************************
 125 
 126   Conjugate gradient solver with partial Cholesky preconditioning
 127   parameter:
 128 
 129            A       compact row oriented storage of lower left of the scaled
 130                    matrix A  
 131            x       solution vector 
 132            b       right hand side
 133            neq     order of A, number of equations 
 134            len     number of non zero entries in Matrix A 
 135            ia      column indices 
 136            iz      row indices (diagonal elements indices)
 137            eps     required accuracy -> residual
 138            niter   maximum number of iterations -> number of iterations
 139            precFlg preconditioning flag 
 140 
 141   The function corresponds to function PACHCG() in H.R. Schwarz: FORTRAN-Pro-
 142   gramme zur Methode der finiten Elemente, p.115, Teubner, 1981
 143 
 144 ********************************************************************** 
 145 */
 146 
 147 void PCG (double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, 
 148           ITG *iz,double *eps, ITG *niter, ITG precFlg,
 149           double *rho, double *r, double *g, double *C, double *z)
 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 }
 249 
 250 
 251 /* **********************************************************************
 252 
 253   Scaling the equation system A x + b = 0
 254 
 255   The equation system is scaled in consideration of keeping the symmetry in
 256   such a way that the diagonal elements of matrix A are 1. This is performed
 257   by the diagonal matrix Ds with the diagonal elements d_i = 1/sqrt(a_ii).
 258   The given equation system Ax+b=0 is transformed into
 259                -1                     - -   -
 260      Ds A Ds  Ds x + Ds b = 0   or    A x + b = 0
 261                          _                                             _ 
 262   with the scaled Matrix A= Ds A Ds  and  the scaled right hand side b= Ds b.
 263   The scaling factor Ds is needed later for backscaling of the solution 
 264   vector
 265      _    -1                       _ 
 266      x = Ds x      resp.    x = Ds x
 267 
 268   parameter:
 269            A       compact row oriented storage of lower left of matrix A
 270            b       right hand side
 271            neq     order of A, number of equations
 272            ia      column indices
 273            iz      row indices (diagonal elements indices)
 274 
 275   The function corresponds to function SCALKO() in H.R. Schwarz: FORTRAN-Pro-
 276   gramme zur Methode der finiten Elemente, p.105, Teubner, 1981
 277 
 278 ********************************************************************** 
 279 */
 280 
 281 void Scaling (double *A, double *b, ITG neq, ITG *ia, ITG *iz, double *d)
 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 }
 305 
 306 /* **********************************************************************
 307 
 308 Computation of matrix vector product z = A p
 309 
 310   parameter:
 311              A       compact row oriented storage of lower left of matrix A
 312              p       vector
 313              neq     order of A, number of equations
 314              ia      column indices
 315              iz      row indices (diagonal elements indices)
 316 
 317   The function corresponds to function APZ() in H.R. Schwarz: FORTRAN-
 318   Programme zur Methode der finiten Elemente, p.111, Teubner, 1981
 319 
 320 ********************************************************************** 
 321 */
 322 
 323 void MatVecProduct (double *A, double *p, ITG neq, ITG *ia, ITG *iz, 
 324                     double *z)
 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 }
 351 
 352 /*--Preconditioning of the equation system Ax+b=0----------------------------------     */
 353 /*---------------------------------------------------------------------------------     */
 354 /*--Partial Cholesky decomposition of scaled matrix A. The off-diagonal matrix   --     */
 355 /*--elements are divided by 1/(1+alpha) until a decomposition of the positive    --     */
 356 /*--definit matrix A exists. C is obtained by ignoring the fill-in during the    --     */
 357 /*--decomposition. In case of successfull decomposition the preconditioning      --     */
 358 /*--matrix C is returned. Otherwise the function is called again with new alpha. --     */
 359 /*--alpha has to be chosen as small as possible, because preconditioning effect  --     */
 360 /*--decreases with increasing alpha.-----------------------------------------------     */
 361 /*---------------------------------------------------------------------------------     */
 362 /*--parameter:                                                                   --     */
 363 /*--           A       compact row oriented storage of lower left of matrix A    --     */
 364 /*--           neq     order of A, number of equations                           --     */
 365 /*--           ia      column indices                                            --     */
 366 /*--           iz      row indices (diagonal elements indices)                   --     */
 367 /*--           alpha   see above                                                 --     */
 368 /*--           precFlg ==3: partial Cholesky partition on full matrix A          --     */
 369 /*---------------------------------------------------------------------------------     */
 370 /*--The function corresponds to function PARTCH() in H.R. Schwarz: FORTRAN-Pro-  --     */
 371 /*--gramme zur Methode der finiten Elemente, p.117, Teubner, 1981                --     */
 372 /*---------------------------------------------------------------------------------     */
 373 void PreConditioning (double *A, ITG neq, ITG len, ITG *ia, ITG *iz, 
 374                                                         double alpha, ITG precFlg,
 375                       double *C, ITG *ier)
 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 }
 425 /*---------------------------------------------------------------------------------     */
 426 
 427 
 428 
 429 /*--Solution of the equation system  M rho = r-------------------------------------     */
 430 /*---------------------------------------------------------------------------------     */
 431 /*--The solution of the equation system M rho = r, with M = C*CT, where C is the --     */
 432 /*--partial Cholesky decomposition of matrix A, represent a preconditioning step.--     */
 433 /*--The equation  system is solved by forward/backward substitution.---------------     */
 434 /*---------------------------------------------------------------------------------     */
 435 /*--parameter:                                                                   --     */
 436 /*--           C       compact row oriented storage of preconditioned matrix A   --     */
 437 /*--           neq     order of A, number of equations                           --     */
 438 /*--           ia      column indices                                            --     */
 439 /*--           iz      row indices (diagonal elements indices)                   --     */
 440 /*--           r       residual vector                                           --     */
 441 /*---------------------------------------------------------------------------------     */
 442 /*--The function corresponds to function CRHOER() in H.R. Schwarz: FORTRAN-Pro-  --     */
 443 /*--gramme zur Methode der finiten Elemente, p.118, Teubner, 1981                --     */
 444 /*---------------------------------------------------------------------------------     */
 445 void Mrhor (double *C, ITG neq, ITG *ia, ITG *iz, double *r, double *rho)
 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 }
 470 /*---------------------------------------------------------------------------------     */
 471 
 472 
 473 
 474 
 475 /*--Calculation of the inner product of two (distributed) vectors------------------     */
 476 /*---------------------------------------------------------------------------------     */
 477 void InnerProduct (double *a, double *b, ITG n, double *Sum)
 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 }
 485 /*---------------------------------------------------------------------------------     */
 486 
 487 
 488 
 489 
 490 /*--(parallel) conjugate gradient solver-------------------------------------------     */
 491 /*---------------------------------------------------------------------------------     */
 492 /*--parameter:                                                                   --     */
 493 /*--           A       compact row oriented storage of lower left of the scaled  --     */
 494 /*--                   matrix A                                                  --     */
 495 /*--           x       solution vector                                           --     */
 496 /*--           b       right hand side                                           --     */
 497 /*--           neq     order of A, number of equations                           --     */
 498 /*--           len     number of non zero entries in Matrix A                    --     */
 499 /*--           ia      column indices                                            --     */
 500 /*--           iz      row indices (diagonal elements indices)                   --     */
 501 /*--           eps     required accuracy -> residual                             --     */
 502 /*--           niter   maximum number of iterations -> number of iterations      --     */
 503 /*---------------------------------------------------------------------------------     */
 504 void CG (double *A, double *x, double *b, ITG neq, ITG len, ITG *ia, ITG *iz,
 505                 double *eps, ITG *niter, double *r, double *p, double *z)
 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 }
 574 /*---------------------------------------------------------------------------------     */

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