# Diff of /trunk/paso/src/SchurComplement.cpp

trunk/paso/src/Solvers/Solver_SchurComplement.c revision 682 by robwdcock, Mon Mar 27 02:43:09 2006 UTC trunk/paso/src/Solver_SchurComplement.c revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1
/* \$Id\$ */
1
2  /*  /*******************************************************
3  ********************************************************************************  *
4  *               Copyright   2006 by ACcESS MNRF                                *  * Copyright (c) 2003-2009 by University of Queensland
5  *                                                                              *  * Earth Systems Science Computational Center (ESSCC)
6  *                 http://www.access.edu.au                                     *  * http://www.uq.edu.au/esscc
7  *           Primary Business: Queensland, Australia                            *  *
11  */  *
12    *******************************************************/
13
14
15  /**************************************************************/  /**************************************************************/
16
# Line 24  Line 25
25
26  /**************************************************************/  /**************************************************************/
27
28  #include "../Paso.h"  #include "Paso.h"
29  #include "../SystemMatrix.h"  #include "SparseMatrix.h"
30  #include "Solver.h"  #include "Solver.h"
31
32  /**************************************************************/  /**************************************************************/
33
34
35
36  void Paso_Solver_updateIncompleteSchurComplement(Paso_SystemMatrix* A_CC,Paso_SystemMatrix *A_CF,double* invA_FF,index_t* A_FF_pivot,Paso_SystemMatrix *A_FC) {  void Paso_Solver_updateIncompleteSchurComplement(Paso_SparseMatrix* A_CC,Paso_SparseMatrix *A_CF,double* invA_FF,index_t* A_FF_pivot,Paso_SparseMatrix *A_FC) {
37
38    index_t iPtr_CC,*index_CC,col_CF,col_FC, *where_p,iPtr_CC_2,i,iPtr_CF,iPtr_FC;    index_t iPtr_CC,*index_CC,col_CF,col_FC, *where_p,iPtr_CC_2,i,iPtr_CF,iPtr_FC;
39    dim_t index_CC_len;    dim_t index_CC_len;
40    bool_t set_A;    bool_t set_A;
41    dim_t n_rows=A_CC->num_rows;    dim_t n_loc_rows=A_CC->numRows;
42    dim_t n_block=A_CC->row_block_size;    dim_t n_block=A_CC->row_block_size;
43    register double A_CF_11,A_CF_21,A_CF_31,A_CF_12,A_CF_22,A_CF_32,A_CF_13,A_CF_23,A_CF_33,    register double A_CF_11,A_CF_21,A_CF_31,A_CF_12,A_CF_22,A_CF_32,A_CF_13,A_CF_23,A_CF_33,
44           invA_FF_11,invA_FF_21,invA_FF_31,invA_FF_12,invA_FF_22,invA_FF_32,invA_FF_13,invA_FF_23,invA_FF_33,           invA_FF_11,invA_FF_21,invA_FF_31,invA_FF_12,invA_FF_22,invA_FF_32,invA_FF_13,invA_FF_23,invA_FF_33,
45           A11,A21,A31,A12,A22,A32,A13,A23,A33,A_FC_11,A_FC_21,A_FC_31,A_FC_12,A_FC_22,A_FC_32,A_FC_13,A_FC_23,A_FC_33;           A11=0,A21=0,A31=0,A12=0,A22=0,A32=0,A13=0,A23=0,A33=0,A_FC_11,A_FC_21,A_FC_31,A_FC_12,A_FC_22,A_FC_32,A_FC_13,A_FC_23,A_FC_33;
46    if (n_block==1) {    if (n_block==1) {
47       #pragma omp parallel for private(i,iPtr_CC,index_CC,index_CC_len,col_CF,set_A,iPtr_CF,iPtr_FC,col_FC,where_p,A11) schedule(static)       #pragma omp parallel for firstprivate(A11) private(i,iPtr_CC,index_CC,index_CC_len,col_CF,set_A,iPtr_CF,iPtr_FC,col_FC,where_p) schedule(static)
48       for (i = 0; i < n_rows;++i) {       for (i = 0; i < n_loc_rows;++i) {
49          iPtr_CC=A_CC->pattern->ptr[i];          iPtr_CC=A_CC->pattern->ptr[i];
50          index_CC=&(A_CC->pattern->index[iPtr_CC]);          index_CC=&(A_CC->pattern->index[iPtr_CC]);
51          index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);          index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);
# Line 54  void Paso_Solver_updateIncompleteSchurCo Line 56  void Paso_Solver_updateIncompleteSchurCo
56               for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) {               for (iPtr_FC = A_FC->pattern->ptr[col_CF]; iPtr_FC < A_FC->pattern->ptr[col_CF + 1]; ++iPtr_FC) {
57                  col_FC=A_FC->pattern->index[iPtr_FC];                  col_FC=A_FC->pattern->index[iPtr_FC];
58                  /* is (i,col_FC) in the shape of A_CC ? */                  /* is (i,col_FC) in the shape of A_CC ? */
59                  where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex);                 where_p=(index_t*)bsearch(&col_FC,index_CC,index_CC_len,sizeof(index_t),Paso_comparIndex);
60                  if (where_p!=NULL) {                  if (where_p!=NULL) {
61                      if (set_A) {                      if (set_A) {
62                         A11=A_CF->val[iPtr_CF]*invA_FF[col_CF];                         A11=A_CF->val[iPtr_CF]*invA_FF[col_CF];
63                         set_A=FALSE;                         set_A=FALSE;
# Line 66  void Paso_Solver_updateIncompleteSchurCo Line 68  void Paso_Solver_updateIncompleteSchurCo
68           } /* end of iPtr_CF loop */           } /* end of iPtr_CF loop */
69        } /* end of irow loop */        } /* end of irow loop */
70     } else if (n_block==2) {     } else if (n_block==2) {
71        #pragma omp parallel for private(i,iPtr_CC,index_CC,index_CC_len,iPtr_CF,col_CF,iPtr_FC,col_FC,where_p,iPtr_CC_2,set_A,A_CF_11,A_CF_21,A_CF_12,A_CF_22,invA_FF_11,invA_FF_21,invA_FF_12,invA_FF_22,A11,A21,A12,A22,A_FC_11,A_FC_21,A_FC_12,A_FC_22) schedule(static)        #pragma omp parallel for firstprivate(A11,A21,A12,A22) private(i,iPtr_CC,index_CC,index_CC_len,iPtr_CF,col_CF,iPtr_FC,col_FC,where_p,iPtr_CC_2,set_A,A_CF_11,A_CF_21,A_CF_12,A_CF_22,invA_FF_11,invA_FF_21,invA_FF_12,invA_FF_22,A_FC_11,A_FC_21,A_FC_12,A_FC_22) schedule(static)
72       for (i = 0; i < n_rows;++i) {       for (i = 0; i < n_loc_rows;++i) {
73          iPtr_CC=A_CC->pattern->ptr[i];          iPtr_CC=A_CC->pattern->ptr[i];
74          index_CC=&(A_CC->pattern->index[iPtr_CC]);          index_CC=&(A_CC->pattern->index[iPtr_CC]);
75          index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);          index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);
# Line 116  void Paso_Solver_updateIncompleteSchurCo Line 118  void Paso_Solver_updateIncompleteSchurCo
118           } /* end of iPtr_CF loop */           } /* end of iPtr_CF loop */
119        } /* end of irow loop */        } /* end of irow loop */
120     } else if (n_block==3) {     } else if (n_block==3) {
121        #pragma omp parallel for private(i,iPtr_CC,index_CC,index_CC_len,iPtr_CF,col_CF,iPtr_FC,col_FC,where_p,iPtr_CC_2,set_A,A_CF_11,A_CF_21,A_CF_31,A_CF_12,A_CF_22,A_CF_32,A_CF_13,A_CF_23,A_CF_33,invA_FF_11,invA_FF_21,invA_FF_31,invA_FF_12,invA_FF_22,invA_FF_32,invA_FF_13,invA_FF_23,invA_FF_33,A11,A21,A31,A12,A22,A32,A13,A23,A33,A_FC_11,A_FC_21,A_FC_31,A_FC_12,A_FC_22,A_FC_32,A_FC_13,A_FC_23,A_FC_33) schedule(static)        #pragma omp parallel for firstprivate(A11,A21,A31,A12,A22,A32,A13,A23,A33) private(i,iPtr_CC,index_CC,index_CC_len,iPtr_CF,col_CF,iPtr_FC,col_FC,where_p,iPtr_CC_2,set_A,A_CF_11,A_CF_21,A_CF_31,A_CF_12,A_CF_22,A_CF_32,A_CF_13,A_CF_23,A_CF_33,invA_FF_11,invA_FF_21,invA_FF_31,invA_FF_12,invA_FF_22,invA_FF_32,invA_FF_13,invA_FF_23,invA_FF_33,A_FC_11,A_FC_21,A_FC_31,A_FC_12,A_FC_22,A_FC_32,A_FC_13,A_FC_23,A_FC_33) schedule(static)
122       for (i = 0; i < n_rows;++i) {       for (i = 0; i < n_loc_rows;++i) {
123          iPtr_CC=A_CC->pattern->ptr[i];          iPtr_CC=A_CC->pattern->ptr[i];
124          index_CC=&(A_CC->pattern->index[iPtr_CC]);          index_CC=&(A_CC->pattern->index[iPtr_CC]);
125          index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);          index_CC_len=(size_t)(A_CC->pattern->ptr[i+1]-A_CC->pattern->ptr[i]);

Legend:
 Removed from v.682 changed lines Added in v.2548