/[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

trunk/paso/src/SparseMatrix_invMain.c revision 3158 by gross, Mon Sep 6 06:09:11 2010 UTC branches/doubleplusgood/paso/src/SparseMatrix_invMain.cpp revision 4261 by jfenwick, Wed Feb 27 06:09:33 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    
17  /**************************************************************/  /************************************************************************************/
18    
19  /* Paso: inverts the main diagonal of a SparseMatrix: */  /* Paso: inverts the main diagonal of a SparseMatrix: */
20    
21  /**************************************************************/  /************************************************************************************/
22    
 /* Copyrights by ACcESS Australia 2010 */  
23  /* Author: Lutz Gross, l.gross@uq.edu.au */  /* Author: Lutz Gross, l.gross@uq.edu.au */
24    
25  /**************************************************************/  /************************************************************************************/
26    
27  #include "Paso.h"  #include "Paso.h"
28  #include "SparseMatrix.h"  #include "SparseMatrix.h"
# Line 30  Line 31 
31  #include "PasoUtil.h"  #include "PasoUtil.h"
32    
33  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) {
34     /*   inv_diag=MEMALLOC( A->numRows * A_p-> block_size,double);     index_t failed=0;
35     pivot=MEMALLOC( A->numRows * A->row_block_size */     register double A11;
36     dim_t n=A_p->numRows;     const dim_t n=A_p->numRows;
37     dim_t n_block=A_p->row_block_size;     const dim_t n_block=A_p->row_block_size;
38     dim_t m_block=A_p->col_block_size;     const dim_t m_block=A_p->col_block_size;
39     /* dim_t block_size=A_p->block_size; */     const dim_t block_size=A_p->block_size;
40     dim_t i;     dim_t i;
41     register index_t iPtr;     register index_t iPtr;
    register double A11,A12,A13,A21,A22,A23,A31,A32,A33,D;  
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.");             failed=1;
58          }          }
59           }           }
60        } else if (n_block==2) {        } else if (n_block==2) {
61       #pragma omp parallel for private(i, iPtr, A11,A12,A21,A22,D) schedule(static)       #pragma omp parallel for private(i, iPtr) schedule(static)
62           for (i = 0; i < n; i++) {           for (i = 0; i < n; i++) {
63             iPtr= main_ptr[i];             iPtr= main_ptr[i];
64             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 {  
           Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");  
            }  
65      }      }
66        } else if (n_block==3) {        } else if (n_block==3) {
67            #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)
68            for (i = 0; i < n; i++) {            for (i = 0; i < n; i++) {
69            iPtr= main_ptr[i];            iPtr= main_ptr[i];
70            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 {  
          Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");  
           }  
71            }            }
72        } else {            } else {    
73       /*       #pragma omp parallel for private(i, iPtr) schedule(static)
      #pragma omp parallel for private(i, INFO) schedule(static)  
74       for (i = 0; i < n; i++) {       for (i = 0; i < n; i++) {
75          iPtr= main_ptr[i];          iPtr= main_ptr[i];
76          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]);
77                    Paso_BlockOps_invM_N(n_block, &inv_diag[i*block_size], &pivot[i*n_block], &failed);
78          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 ) {  
            Paso_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");  
           }  
      */  
      Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_invMain: Right now there is support block size less than 4 only");  
79        }        }
80     }     }
81       if (failed > 0) {
82          Esys_setError(ZERO_DIVISION_ERROR, "Paso_SparseMatrix_invMain: non-regular main diagonal block.");
83       }
84  }  }
85  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) {
86        
# Line 131  void Paso_SparseMatrix_applyBlockMatrix( Line 89  void Paso_SparseMatrix_applyBlockMatrix(
89     dim_t n=A_p->numRows;     dim_t n=A_p->numRows;
90     dim_t n_block=A_p->row_block_size;     dim_t n_block=A_p->row_block_size;
91     Paso_Copy(n_block*n, x,b);     Paso_Copy(n_block*n, x,b);
92     Paso_BlockOps_allMV(n_block,n,block_diag,pivot,x);     Paso_BlockOps_solveAll(n_block,n,block_diag,pivot,x);
93  }  }
94    

Legend:
Removed from v.3158  
changed lines
  Added in v.4261

  ViewVC Help
Powered by ViewVC 1.1.26