/[escript]/trunk/paso/src/Solver_jacobi.c
ViewVC logotype

Diff of /trunk/paso/src/Solver_jacobi.c

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

revision 3093 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 3094 by gross, Fri Aug 13 08:38:06 2010 UTC
# Line 37  void Paso_Solver_Jacobi_free(Paso_Solver Line 37  void Paso_Solver_Jacobi_free(Paso_Solver
37       MEMFREE(in);       MEMFREE(in);
38    }    }
39  }  }
40    Paso_Solver_Jacobi* Paso_Solver_getJacobi(Paso_SystemMatrix * A_p)
41    {
42       Paso_Solver_Jacobi* jac=Paso_Solver_getLocalJacobi(A_p->mainBlock);
43       if (Paso_MPIInfo_noError(A_p->mpi_info)) {
44             return jac;
45       } else {
46             return NULL;
47       }
48    }
49    
50  /**************************************************************/  /**************************************************************/
51    
52  /* Jacobi precondioner set up */  /* Jacobi precondioner set up */
53    
54  Paso_Solver_Jacobi* Paso_Solver_getJacobi(Paso_SparseMatrix * A_p) {  Paso_Solver_Jacobi* Paso_Solver_getLocalJacobi(Paso_SparseMatrix * A_p) {
55    Paso_Solver_Jacobi* out=NULL;    Paso_Solver_Jacobi* out=NULL;
56    dim_t n = A_p->numCols;    dim_t n = A_p->numCols;
57    dim_t n_block=A_p->row_block_size;    dim_t n_block=A_p->row_block_size;
58    dim_t block_size=A_p->block_size;    dim_t block_size=A_p->block_size;
59    dim_t i;    
   index_t iPtr;  
   double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;  
60    /* check matrix is square */    /* check matrix is square */
61    if (A_p->col_block_size !=A_p->row_block_size) {    if (A_p->col_block_size !=A_p->row_block_size) {
62      Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Jacobi preconditioner square block size.");      Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Jacobi preconditioner square block size.");
63      return NULL;      return NULL;
64    }    }
65    /* check matrix is square */    
   if (n_block>3) {  
     Paso_setError(TYPE_ERROR, "Paso_Solver_getJacobi: Right now the Jacobi preconditioner supports block size less than 4 only");  
     return NULL;  
   }  
   /* allocate vector to hold main diagonal entries: */  
66    out=MEMALLOC(1,Paso_Solver_Jacobi);    out=MEMALLOC(1,Paso_Solver_Jacobi);
67    if (! Paso_checkPtr(out)) {    if (! Paso_checkPtr(out)) {
68        /* allocate vector to hold main diagonal entries: */       out->values = MEMALLOC( ((size_t) n) * ((size_t) block_size),double);
69        out->n_block=n_block;       out->pivot = MEMALLOC(  ((size_t) n) * ((size_t) n_block), index_t);
70        out->n=n;  
71        out->values = MEMALLOC( ((size_t) n) * ((size_t) block_size),double);      
72        out->pivot = NULL; /* later use */       if (! ( Paso_checkPtr(out->values) || Paso_checkPtr(out->pivot) )) {
73        if (! (Paso_checkPtr(out->values))) {           Paso_SparseMatrix_invMain(A_p, out->values, out->pivot );
74          if (n_block==1) {       }
75             #pragma omp parallel for private(i, iPtr) schedule(static)      
            for (i = 0; i < A_p->pattern->numOutput; i++) {  
               out->values[i]=1.;  
               /* find main diagonal */  
               for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {  
                   if (A_p->pattern->index[iPtr]==i) {  
                       if (ABS(A_p->val[iPtr])>0.) out->values[i]=1./A_p->val[iPtr];  
                       break;  
                   }  
               }  
            }  
         } else if (n_block==2) {  
            #pragma omp parallel for private(i, iPtr, A11,A12,A21,A22,D) schedule(static)  
            for (i = 0; i < A_p->pattern->numOutput; i++) {  
               out->values[i*4+0]= 1.;  
               out->values[i*4+1]= 0.;  
               out->values[i*4+2]= 0.;  
               out->values[i*4+3]= 1.;  
               /* find main diagonal */  
               for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {  
                   if (A_p->pattern->index[iPtr]==i) {  
                      A11=A_p->val[iPtr*4];  
                      A12=A_p->val[iPtr*4+2];  
                      A21=A_p->val[iPtr*4+1];  
                      A22=A_p->val[iPtr*4+3];  
                      D = A11*A22-A12*A21;  
                      if (ABS(D) > 0 ){  
                         D=1./D;  
                         out->values[i*4]= A22*D;  
                         out->values[i*4+1]=-A21*D;  
                         out->values[i*4+2]=-A12*D;  
                         out->values[i*4+3]= A11*D;  
                      }  
                      break;  
                   }  
               }  
            }  
         } else if (n_block==3) {  
            #pragma omp parallel for private(i, iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D) schedule(static)  
            for (i = 0; i < A_p->pattern->numOutput; i++) {  
               out->values[i*9  ]=1.;  
               out->values[i*9+1]=0.;  
               out->values[i*9+2]=0.;  
               out->values[i*9+3]=0.;  
               out->values[i*9+4]=1.;  
               out->values[i*9+5]=0.;  
               out->values[i*9+6]=0.;  
               out->values[i*9+7]=0.;  
               out->values[i*9+8]=1.;  
               /* find main diagonal */  
               for (iPtr = A_p->pattern->ptr[i]; iPtr < A_p->pattern->ptr[i + 1]; iPtr++) {  
                   if (A_p->pattern->index[iPtr]==i) {  
                       A11=A_p->val[iPtr*9  ];  
                       A21=A_p->val[iPtr*9+1];  
                       A31=A_p->val[iPtr*9+2];  
                       A12=A_p->val[iPtr*9+3];  
                       A22=A_p->val[iPtr*9+4];  
                       A32=A_p->val[iPtr*9+5];  
                       A13=A_p->val[iPtr*9+6];  
                       A23=A_p->val[iPtr*9+7];  
                       A33=A_p->val[iPtr*9+8];  
                       D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);  
                       if (ABS(D) > 0 ){  
                          D=1./D;  
                          out->values[i*9  ]=(A22*A33-A23*A32)*D;  
                          out->values[i*9+1]=(A31*A23-A21*A33)*D;  
                          out->values[i*9+2]=(A21*A32-A31*A22)*D;  
                          out->values[i*9+3]=(A13*A32-A12*A33)*D;  
                          out->values[i*9+4]=(A11*A33-A31*A13)*D;  
                          out->values[i*9+5]=(A12*A31-A11*A32)*D;  
                          out->values[i*9+6]=(A12*A23-A13*A22)*D;  
                          out->values[i*9+7]=(A13*A21-A11*A23)*D;  
                          out->values[i*9+8]=(A11*A22-A12*A21)*D;  
                       }  
                       break;  
                   }  
               }  
            }  
         }  
       }  
76    }    }
77    if (Paso_noError()) {    if (Paso_noError()) {
78       return out;       return out;
# Line 165  Paso_Solver_Jacobi* Paso_Solver_getJacob Line 88  Paso_Solver_Jacobi* Paso_Solver_getJacob
88  /* should be called within a parallel region                                              */  /* should be called within a parallel region                                              */
89  /* 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 */
90    
91  void Paso_Solver_solveJacobi(Paso_Solver_Jacobi * prec, double * x, double * b) {  void Paso_Solver_solveJacobi(Paso_SystemMatrix * A_p, Paso_Solver_Jacobi * prec, double * x, double * b) {
92       Paso_Solver_applyBlockDiagonalMatrix(prec->n_block,prec->n,prec->values,prec->pivot,x,b);     Paso_Solver_solveLocalJacobi( A_p->mainBlock, prec, x,  b) ;
93       return;
94    }
95    
96    void Paso_Solver_solveLocalJacobi(Paso_SparseMatrix * A_p, Paso_Solver_Jacobi * prec, double * x, double * b) {
97         dim_t n = A_p->numCols;
98         dim_t n_block=A_p->row_block_size;
99         Paso_Solver_applyBlockDiagonalMatrix(n_block,n,prec->values,prec->pivot,x,b);
100       return;       return;
101  }  }
102    

Legend:
Removed from v.3093  
changed lines
  Added in v.3094

  ViewVC Help
Powered by ViewVC 1.1.26