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

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

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

revision 3258 by gross, Mon Sep 6 06:09:11 2010 UTC revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 42  void Paso_SparseMatrix_invMain(Paso_Spar Line 42  void Paso_SparseMatrix_invMain(Paso_Spar
42     index_t* main_ptr=Paso_Pattern_borrowMainDiagonalPointer(A_p->pattern);     index_t* main_ptr=Paso_Pattern_borrowMainDiagonalPointer(A_p->pattern);
43     /* check matrix is square */     /* check matrix is square */
44     if (m_block != n_block) {     if (m_block != n_block) {
45        Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_invMain: square block size expected.");        Esys_setError(TYPE_ERROR, "Paso_SparseMatrix_invMain: square block size expected.");
46     }     }
47     if (Paso_noError()) {     if (Esys_noError()) {
48                    
49        if (n_block==1) {        if (n_block==1) {
50            #pragma omp parallel for private(i, iPtr, A11) schedule(static)            #pragma omp parallel for private(i, iPtr, A11) schedule(static)
# Line 54  void Paso_SparseMatrix_invMain(Paso_Spar Line 54  void Paso_SparseMatrix_invMain(Paso_Spar
54          if ( ABS(A11) > 0.) {          if ( ABS(A11) > 0.) {
55              inv_diag[i]=1./A11;              inv_diag[i]=1./A11;
56          } else {          } else {
57             Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");             Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");
58          }          }
59           }           }
60        } else if (n_block==2) {        } else if (n_block==2) {
# Line 73  void Paso_SparseMatrix_invMain(Paso_Spar Line 73  void Paso_SparseMatrix_invMain(Paso_Spar
73           inv_diag[i*4+2]=-A12*D;           inv_diag[i*4+2]=-A12*D;
74           inv_diag[i*4+3]= A11*D;           inv_diag[i*4+3]= A11*D;
75             } else {             } else {
76            Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");            Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");
77             }             }
78      }      }
79        } else if (n_block==3) {        } else if (n_block==3) {
# Line 102  void Paso_SparseMatrix_invMain(Paso_Spar Line 102  void Paso_SparseMatrix_invMain(Paso_Spar
102           inv_diag[i*9+7]=(A13*A21-A11*A23)*D;           inv_diag[i*9+7]=(A13*A21-A11*A23)*D;
103           inv_diag[i*9+8]=(A11*A22-A12*A21)*D;           inv_diag[i*9+8]=(A11*A22-A12*A21)*D;
104            } else {            } else {
105           Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");           Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");
106            }            }
107            }            }
108        } else {            } else {    
# Line 117  void Paso_SparseMatrix_invMain(Paso_Spar Line 117  void Paso_SparseMatrix_invMain(Paso_Spar
117                    
118          dgetrf_(n_block, m_block, &(diag[i*block_size], n_block, pivot[i*n_block], info );          dgetrf_(n_block, m_block, &(diag[i*block_size], n_block, pivot[i*n_block], info );
119          if (INFO >0 ) {          if (INFO >0 ) {
120             Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");             Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");
121            }            }
122       */       */
123       Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_invMain: Right now there is support block size less than 4 only");       Esys_setError(TYPE_ERROR, "Paso_SparseMatrix_invMain: Right now there is support block size less than 4 only");
124        }        }
125     }     }
126  }  }

Legend:
Removed from v.3258  
changed lines
  Added in v.3259

  ViewVC Help
Powered by ViewVC 1.1.26