root/src/preiter.c

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

DEFINITIONS

This source file includes following definitions.
  1. preiter

   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){
  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 }

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