/[escript]/trunk/esys2/finley/src/finleyC/Solvers/Solver_jacobi.c
ViewVC logotype

Diff of /trunk/esys2/finley/src/finleyC/Solvers/Solver_jacobi.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 114 by jgs, Wed Dec 15 07:08:39 2004 UTC revision 115 by jgs, Fri Mar 4 07:12:47 2005 UTC
# Line 6  Line 6 
6    
7  /**************************************************************/  /**************************************************************/
8    
9  /* Copyrights by ACcESS Australia 2003/04 */  /* Copyrights by ACcESS Australia 2003, 2004, 2005 */
10  /* Author: gross@access.edu.au */  /* Author: gross@access.edu.au */
11    
12  /**************************************************************/  /**************************************************************/
13    
 #include "Finley.h"  
 #include "System.h"  
14  #include "Solver.h"  #include "Solver.h"
 #include "Util.h"  
 #include "Common.h"  
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /* inverse preconditioner setup */  /* frees the Jacobi preconditioner */
19    
20  void Finley_Solver_setJacobi(Finley_SystemMatrix * A_p) {  void Finley_Solver_Jacobi_free(Finley_Solver_Jacobi * in) {
21  #if ITERATIVE_SOLVER == NO_LIB    if (in!=NULL) {
22    Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative);       MEMFREE(in->values);
23         MEMFREE(in->pivot);
24         MEMFREE(in);
25      }
26    }
27    
28    /**************************************************************/
29    
30    /* Jacobi precondioner set up */
31    
32    Finley_Solver_Jacobi* Finley_Solver_getJacobi(Finley_SystemMatrix * A_p) {
33      Finley_Solver_Jacobi* out=NULL;
34    int n = A_p->num_cols;    int n = A_p->num_cols;
35    int n_block=A_p->row_block_size;    int n_block=A_p->row_block_size;
36    int block_size=A_p->block_size;    int block_size=A_p->block_size;
37    int i,iPtr,info;    int i,iPtr;
38      double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;
39    /* check matrix is square */    /* check matrix is square */
40    if (A_p->col_block_size !=A_p->row_block_size) {    if (A_p->col_block_size !=A_p->row_block_size) {
41      Finley_ErrorCode = TYPE_ERROR;      Finley_ErrorCode = TYPE_ERROR;
42      sprintf(Finley_ErrorMsg, "Jacobi preconditioner square block size.");      sprintf(Finley_ErrorMsg, "Jacobi preconditioner square block size.");
43      return;      return NULL;
44    }    }
45    /* check matrix is square */    /* check matrix is square */
46    if (n_block>3) {    if (n_block>3) {
47      Finley_ErrorCode = TYPE_ERROR;      Finley_ErrorCode = TYPE_ERROR;
48      sprintf(Finley_ErrorMsg, "Right now the Jacobi preconditioner supports block size less than 4 only");      sprintf(Finley_ErrorMsg, "Right now the Jacobi preconditioner supports block size less than 4 only");
49      return;      return NULL;
50    }    }
51    /* allocate vector to hold main diagonal entries: */    /* allocate vector to hold main diagonal entries: */
52    prec->values = MEMALLOC( ((size_t) n) * ((size_t) block_size),double);    out=MEMALLOC(1,Finley_Solver_Jacobi);
53    prec->pivot = MEMALLOC(  ((size_t) n) * ((size_t) n_block),int);    if (! Finley_checkPtr(out)) {
54    if (! (Finley_checkPtr(prec->values) || Finley_checkPtr(prec->pivot)) ) {        /* allocate vector to hold main diagonal entries: */
55          out->n_block=n_block;
56       if (n_block==1) {        out->n=n;
57          #pragma omp parallel for private(i, iPtr) schedule(static)        out->values = MEMALLOC( ((size_t) n) * ((size_t) block_size),double);
58          for (i = 0; i < A_p->pattern->n_ptr; i++) {        out->pivot = NULL; /* later use */
59             /* find main diagonal */        if (! (Finley_checkPtr(out->values))) {
60             for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {          if (n_block==1) {
61                 if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) {             #pragma omp parallel for private(i, iPtr) schedule(static)
62                     if (ABS(A_p->val[iPtr])>0.) {             for (i = 0; i < A_p->pattern->n_ptr; i++) {
63                        prec->values[i]=1./A_p->val[iPtr];                out->values[i]=1.;
64                     } else {                /* find main diagonal */
65                        prec->values[i]=1.;                for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {
66                     }                    if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) {
67                     break;                        if (ABS(A_p->val[iPtr])>0.) out->values[i]=1./A_p->val[iPtr];
68                 }                        break;
69                      }
70                  }
71             }             }
72          }          } else if (n_block==2) {
73       } else {             #pragma omp parallel for private(i, iPtr, A11,A12,A21,A22,D) schedule(static)
74          #pragma omp parallel for private(i, iPtr,info) schedule(static)             for (i = 0; i < A_p->pattern->n_ptr; i++) {
75          for (i = 0; i < A_p->pattern->n_ptr; i++) {                out->values[i*4+0]= 1.;
76             /* find main diagonal */                out->values[i*4+1]= 0.;
77             for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {                out->values[i*4+2]= 0.;
78                 if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) {                out->values[i*4+3]= 1.;
79                     info=Finley_Util_SmallMatLU(n_block,&A_p->val[iPtr*block_size],&prec->values[i*block_size],&prec->pivot[i*n_block]);                /* find main diagonal */
80                     if (info>0) {                for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {
81                          Finley_ErrorCode = ZERO_DIVISION_ERROR;                    if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) {
82                          sprintf(Finley_ErrorMsg, "non-regular main diagonal block in row %d",i);                       A11=A_p->val[iPtr*4];
83                     }                       A12=A_p->val[iPtr*4+2];
84                     break;                       A21=A_p->val[iPtr*4+1];
85                 }                       A22=A_p->val[iPtr*4+3];
86                         D = A11*A22-A12*A21;
87                         if (ABS(D) > 0 ){
88                            D=1./D;
89                            out->values[i*4]= A22*D;
90                            out->values[i*4+1]=-A21*D;
91                            out->values[i*4+2]=-A12*D;
92                            out->values[i*4+3]= A11*D;
93                         }
94                         break;
95                      }
96                  }
97               }
98            } else if (n_block==3) {
99               #pragma omp parallel for private(i, iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D) schedule(static)
100               for (i = 0; i < A_p->pattern->n_ptr; i++) {
101                  out->values[i*9  ]=1.;
102                  out->values[i*9+1]=0.;
103                  out->values[i*9+2]=0.;
104                  out->values[i*9+3]=0.;
105                  out->values[i*9+4]=1.;
106                  out->values[i*9+5]=0.;
107                  out->values[i*9+6]=0.;
108                  out->values[i*9+7]=0.;
109                  out->values[i*9+8]=1.;
110                  /* find main diagonal */
111                  for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {
112                      if (A_p->pattern->index[iPtr]==i+INDEX_OFFSET) {
113                          A11=A_p->val[iPtr*9  ];
114                          A21=A_p->val[iPtr*9+1];
115                          A31=A_p->val[iPtr*9+2];
116                          A12=A_p->val[iPtr*9+3];
117                          A22=A_p->val[iPtr*9+4];
118                          A32=A_p->val[iPtr*9+5];
119                          A13=A_p->val[iPtr*9+6];
120                          A23=A_p->val[iPtr*9+7];
121                          A33=A_p->val[iPtr*9+8];
122                          D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
123                          if (ABS(D) > 0 ){
124                             D=1./D;
125                             out->values[i*9  ]=(A22*A33-A23*A32)*D;
126                             out->values[i*9+1]=(A31*A23-A21*A33)*D;
127                             out->values[i*9+2]=(A21*A32-A31*A22)*D;
128                             out->values[i*9+3]=(A13*A32-A12*A33)*D;
129                             out->values[i*9+4]=(A11*A33-A31*A13)*D;
130                             out->values[i*9+5]=(A12*A31-A11*A32)*D;
131                             out->values[i*9+6]=(A12*A23-A13*A22)*D;
132                             out->values[i*9+7]=(A13*A21-A11*A23)*D;
133                             out->values[i*9+8]=(A11*A22-A12*A21)*D;
134                          }
135                          break;
136                      }
137                  }
138             }             }
139          }          }
140       }        }
141      }
142      if (Finley_ErrorCode == NO_ERROR) {
143         return out;
144      } else {
145         Finley_Solver_Jacobi_free(out);
146         return NULL;
147    }    }
 #endif  
148  }  }
   
149  /**************************************************************/  /**************************************************************/
150    
151  /* inverse preconditioner solution using raw pointers */  /* applies Jacobi preconditioner */
152    
153  /* should be called within a parallel region                                              */  /* should be called within a parallel region                                              */
154  /* barrier synconization should be performed to make sure that the input vector available */  /* barrier synconization should be performed to make sure that the input vector available */
155    
156  void Finley_Solver_solveJacobi(Finley_SystemMatrix * A_p, double * x, double * b) {  void Finley_Solver_solveJacobi(Finley_Solver_Jacobi * prec, double * x, double * b) {
157  #if ITERATIVE_SOLVER == NO_LIB       Finley_Solver_applyBlockDiagonalMatrix(prec->n_block,prec->n,prec->values,prec->pivot,x,b);
158    Finley_Solver_Preconditioner* prec=(Finley_Solver_Preconditioner*) (A_p->iterative);       return;
   int n_block=A_p->row_block_size;  
   int block_size=A_p->block_size;  
   int i;  
   
   if (n_block==1) {  
      #pragma omp for private(i) schedule(static)  
      for (i = 0; i < A_p->pattern->n_ptr; i++) {  
         x[i] = prec->values[i] * b[i];  
      }  
   } else {  
      #pragma omp for private(i) schedule(static)  
      for (i = 0; i < A_p->pattern->n_ptr; i++) {  
         Finley_Util_SmallMatForwardBackwardSolve(n_block,1,&prec->values[i*block_size],&prec->pivot[i*n_block],&x[i*n_block],&b[i*n_block]);  
      }  
   }  
   return;  
 #endif  
159  }  }

Legend:
Removed from v.114  
changed lines
  Added in v.115

  ViewVC Help
Powered by ViewVC 1.1.26