/[escript]/branches/doubleplusgood/paso/src/SparseMatrix_invMain.cpp
ViewVC logotype

Diff of /branches/doubleplusgood/paso/src/SparseMatrix_invMain.cpp

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

revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC revision 3402 by gross, Tue Dec 7 07:36:12 2010 UTC
# Line 30  Line 30 
30  #include "PasoUtil.h"  #include "PasoUtil.h"
31    
32  void Paso_SparseMatrix_invMain(Paso_SparseMatrix * A_p, double* inv_diag, int* pivot) {  void Paso_SparseMatrix_invMain(Paso_SparseMatrix * A_p, double* inv_diag, int* pivot) {
33     /*   inv_diag=MEMALLOC( A->numRows * A_p-> block_size,double);     index_t failed=0;
34     pivot=MEMALLOC( A->numRows * A->row_block_size */     register double A11;
35     dim_t n=A_p->numRows;     const dim_t n=A_p->numRows;
36     dim_t n_block=A_p->row_block_size;     const dim_t n_block=A_p->row_block_size;
37     dim_t m_block=A_p->col_block_size;     const dim_t m_block=A_p->col_block_size;
38     /* dim_t block_size=A_p->block_size; */     const dim_t block_size=A_p->block_size;
39     dim_t i;     dim_t i;
40     register index_t iPtr;     register index_t iPtr;
    register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;  
41     index_t* main_ptr=Paso_Pattern_borrowMainDiagonalPointer(A_p->pattern);     index_t* main_ptr=Paso_Pattern_borrowMainDiagonalPointer(A_p->pattern);
42     /* check matrix is square */     /* check matrix is square */
43     if (m_block != n_block) {     if (m_block != n_block) {
# Line 54  void Paso_SparseMatrix_invMain(Paso_Spar Line 53  void Paso_SparseMatrix_invMain(Paso_Spar
53          if ( ABS(A11) > 0.) {          if ( ABS(A11) > 0.) {
54              inv_diag[i]=1./A11;              inv_diag[i]=1./A11;
55          } else {          } else {
56             Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");             failed=1;
57          }          }
58           }           }
59        } else if (n_block==2) {        } else if (n_block==2) {
60       #pragma omp parallel for private(i, iPtr, A11,A12,A21,A22,D) schedule(static)       #pragma omp parallel for private(i, iPtr) schedule(static)
61           for (i = 0; i < n; i++) {           for (i = 0; i < n; i++) {
62             iPtr= main_ptr[i];             iPtr= main_ptr[i];
63             A11=A_p->val[iPtr*4];             Paso_BlockOps_invM_2(&inv_diag[i*4], &A_p->val[iPtr*4], &failed);
            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;  
          inv_diag[i*4]= A22*D;  
          inv_diag[i*4+1]=-A21*D;  
          inv_diag[i*4+2]=-A12*D;  
          inv_diag[i*4+3]= A11*D;  
            } else {  
           Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");  
            }  
64      }      }
65        } else if (n_block==3) {        } else if (n_block==3) {
66            #pragma omp parallel for private(i, iPtr,A11,A12,A13,A21,A22,A23,A31,A32,A33,D) schedule(static)            #pragma omp parallel for private(i, iPtr) schedule(static)
67            for (i = 0; i < n; i++) {            for (i = 0; i < n; i++) {
68            iPtr= main_ptr[i];            iPtr= main_ptr[i];
69            A11=A_p->val[iPtr*9  ];            Paso_BlockOps_invM_3(&inv_diag[i*9], &A_p->val[iPtr*9], &failed);
           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;  
          inv_diag[i*9  ]=(A22*A33-A23*A32)*D;  
          inv_diag[i*9+1]=(A31*A23-A21*A33)*D;  
          inv_diag[i*9+2]=(A21*A32-A31*A22)*D;  
          inv_diag[i*9+3]=(A13*A32-A12*A33)*D;  
          inv_diag[i*9+4]=(A11*A33-A31*A13)*D;  
          inv_diag[i*9+5]=(A12*A31-A11*A32)*D;  
          inv_diag[i*9+6]=(A12*A23-A13*A22)*D;  
          inv_diag[i*9+7]=(A13*A21-A11*A23)*D;  
          inv_diag[i*9+8]=(A11*A22-A12*A21)*D;  
           } else {  
          Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");  
           }  
70            }            }
71        } else {            } else {    
72       /*       #pragma omp parallel for private(i, iPtr) schedule(static)
      #pragma omp parallel for private(i, INFO) schedule(static)  
73       for (i = 0; i < n; i++) {       for (i = 0; i < n; i++) {
74          iPtr= main_ptr[i];          iPtr= main_ptr[i];
75          memcpy((void*)&(diag[i*block_size], (void*)(&A_p->val[iPtr*block_size]), sizeof(double)*(size_t)block_size);          Paso_BlockOps_Cpy_N(block_size, &inv_diag[i*block_size], &A_p->val[iPtr*block_size]);
76                    Paso_BlockOps_invM_N(n_block, &inv_diag[i*block_size], &pivot[i*n_block], &failed);
77          int dgetrf_(integer *m, integer *n, doublereal *a, integer *       }
         lda, integer *ipiv, integer *info);  
           
         dgetrf_(n_block, m_block, &(diag[i*block_size], n_block, pivot[i*n_block], info );  
         if (INFO >0 ) {  
            Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");  
           }  
      */  
      Esys_setError(TYPE_ERROR, "Paso_SparseMatrix_invMain: Right now there is support block size less than 4 only");  
78        }        }
79     }     }
80       if (failed > 0) {
81          Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");
82       }
83  }  }
84  void Paso_SparseMatrix_applyBlockMatrix(Paso_SparseMatrix * A_p, double* block_diag, int* pivot, double*x, double *b) {  void Paso_SparseMatrix_applyBlockMatrix(Paso_SparseMatrix * A_p, double* block_diag, int* pivot, double*x, double *b) {
85        
# Line 131  void Paso_SparseMatrix_applyBlockMatrix( Line 88  void Paso_SparseMatrix_applyBlockMatrix(
88     dim_t n=A_p->numRows;     dim_t n=A_p->numRows;
89     dim_t n_block=A_p->row_block_size;     dim_t n_block=A_p->row_block_size;
90     Paso_Copy(n_block*n, x,b);     Paso_Copy(n_block*n, x,b);
91     Paso_BlockOps_allMV(n_block,n,block_diag,pivot,x);     Paso_BlockOps_solveAll(n_block,n,block_diag,pivot,x);
92  }  }
93    

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

  ViewVC Help
Powered by ViewVC 1.1.26