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

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

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

revision 700 by gross, Thu Apr 6 00:13:40 2006 UTC revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC
# Line 1  Line 1 
 /* $Id$ */  
1    
2    /* $Id$ */
3    
4  /*  /*******************************************************
5  ********************************************************************************   *
6  *               Copyright   2006 by ACcESS MNRF                                *   *           Copyright 2003-2007 by ACceSS MNRF
7  *                                                                              *   *       Copyright 2007 by University of Queensland
8  *                 http://www.access.edu.au                                     *   *
9  *           Primary Business: Queensland, Australia                            *   *                http://esscc.uq.edu.au
10  *     Licensed under the Open Software License version 3.0             *   *        Primary Business: Queensland, Australia
11  *        http://www.opensource.org/licenses/osl-3.0.php                        *   *  Licensed under the Open Software License version 3.0
12  ********************************************************************************   *     http://www.opensource.org/licenses/osl-3.0.php
13  */   *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
# Line 18  Line 19 
19    
20  /**************************************************************/  /**************************************************************/
21    
22  /* Copyrights by ACcESS Australia 2003,2004,2005              */  /* Copyrights by ACcESS Australia 2004,2004,2005, 2006, 2007  */
23  /* Author: gross@access.edu.au                                */  /* Author: gross@access.edu.au                                */
24    
25  /**************************************************************/  /**************************************************************/
# Line 36  void Paso_Solver_RILU_free(Paso_Solver_R Line 37  void Paso_Solver_RILU_free(Paso_Solver_R
37          Paso_Solver_RILU_free(in->RILU_of_Schur);          Paso_Solver_RILU_free(in->RILU_of_Schur);
38          MEMFREE(in->inv_A_FF);          MEMFREE(in->inv_A_FF);
39          MEMFREE(in->A_FF_pivot);          MEMFREE(in->A_FF_pivot);
40          Paso_SystemMatrix_dealloc(in->A_FC);          Paso_SparseMatrix_free(in->A_FC);
41          Paso_SystemMatrix_dealloc(in->A_CF);          Paso_SparseMatrix_free(in->A_CF);
42          MEMFREE(in->rows_in_F);          MEMFREE(in->rows_in_F);
43          MEMFREE(in->rows_in_C);          MEMFREE(in->rows_in_C);
44          MEMFREE(in->mask_F);          MEMFREE(in->mask_F);
# Line 69  to Line 70  to
70     then RILU is applied to S again until S becomes empty     then RILU is applied to S again until S becomes empty
71    
72  */  */
73  Paso_Solver_RILU* Paso_Solver_getRILU(Paso_SystemMatrix * A_p,bool_t verbose) {  Paso_Solver_RILU* Paso_Solver_getRILU(Paso_SparseMatrix *A_p,bool_t verbose) {
74    Paso_Solver_RILU* out=NULL;    Paso_Solver_RILU* out=NULL;
75    dim_t n=A_p->num_rows;    dim_t n=A_p->numRows;
76    dim_t n_block=A_p->row_block_size;    dim_t n_block=A_p->row_block_size;
77    index_t* mis_marker=NULL;      index_t* mis_marker=NULL;  
78    index_t* counter=NULL;      index_t* counter=NULL;  
79    index_t iPtr,*index, *where_p;    index_t iPtr,*index, *where_p;
80    dim_t i,k;    dim_t i,k;
81    Paso_SystemMatrix * schur=NULL;    Paso_SparseMatrix * schur=NULL;
82    double A11,A12,A13,A21,A22,A23,A31,A32,A33,D,time0,time1,time2;    double A11,A12,A13,A21,A22,A23,A31,A32,A33,D,time0,time1,time2;
83        
84    
# Line 104  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 105  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
105       time0=Paso_timer();       time0=Paso_timer();
106       #pragma omp parallel for private(i) schedule(static)       #pragma omp parallel for private(i) schedule(static)
107       for (i=0;i<n;++i) mis_marker[i]=-1;       for (i=0;i<n;++i) mis_marker[i]=-1;
108       Paso_SystemMatrixPattern_mis(A_p->pattern,mis_marker);       Paso_Pattern_mis(A_p->pattern,mis_marker);
109       time2=Paso_timer()-time0;       time2=Paso_timer()-time0;
110       if (Paso_noError()) {       if (Paso_noError()) {
111          #pragma omp parallel for private(i) schedule(static)          #pragma omp parallel for private(i) schedule(static)
# Line 223  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 224  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
224                            }                            }
225                        } /* end parallel region */                        } /* end parallel region */
226                        /* get A_CF block: */                        /* get A_CF block: */
227                        out->A_CF=Paso_SystemMatrix_getSubmatrix(A_p,out->n_C,out->rows_in_C,out->mask_F);                        out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F);
228                        if (Paso_noError()) {                        if (Paso_noError()) {
229                           /* get A_FC block: */                           /* get A_FC block: */
230                           out->A_FC=Paso_SystemMatrix_getSubmatrix(A_p,out->n_F,out->rows_in_F,out->mask_C);                           out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);
231                           /* get A_FF block: */                           /* get A_FF block: */
232                           if (Paso_noError()) {                           if (Paso_noError()) {
233                              schur=Paso_SystemMatrix_getSubmatrix(A_p,out->n_C,out->rows_in_C,out->mask_C);                              schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);
234                              time0=Paso_timer()-time0;                              time0=Paso_timer()-time0;
235                              if (Paso_noError()) {                              if (Paso_noError()) {
236                                  time1=Paso_timer();                                  time1=Paso_timer();
# Line 237  Paso_Solver_RILU* Paso_Solver_getRILU(Pa Line 238  Paso_Solver_RILU* Paso_Solver_getRILU(Pa
238                                  Paso_Solver_updateIncompleteSchurComplement(schur,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);                                  Paso_Solver_updateIncompleteSchurComplement(schur,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);
239                                  time1=Paso_timer()-time1;                                  time1=Paso_timer()-time1;
240                                  out->RILU_of_Schur=Paso_Solver_getRILU(schur,verbose);                                  out->RILU_of_Schur=Paso_Solver_getRILU(schur,verbose);
241                                  Paso_SystemMatrix_dealloc(schur);                                  Paso_SparseMatrix_free(schur);
242                              }                              }
243                              /* allocate work arrays for RILU application */                              /* allocate work arrays for RILU application */
244                              if (Paso_noError()) {                              if (Paso_noError()) {
# Line 340  void Paso_Solver_solveRILU(Paso_Solver_R Line 341  void Paso_Solver_solveRILU(Paso_Solver_R
341          /* x_F=invA_FF*b_F  */          /* x_F=invA_FF*b_F  */
342          Paso_Solver_applyBlockDiagonalMatrix(n_block,rilu->n_F,rilu->inv_A_FF,rilu->A_FF_pivot,rilu->x_F,rilu->b_F);          Paso_Solver_applyBlockDiagonalMatrix(n_block,rilu->n_F,rilu->inv_A_FF,rilu->A_FF_pivot,rilu->x_F,rilu->b_F);
343          /* b_C=b_C-A_CF*x_F */          /* b_C=b_C-A_CF*x_F */
344          Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(-1.,rilu->A_CF,rilu->x_F,1.,rilu->b_C);          Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,rilu->A_CF,rilu->x_F,1.,rilu->b_C);
345          /* x_C=RILU(b_C)     */          /* x_C=RILU(b_C)     */
346          Paso_Solver_solveRILU(rilu->RILU_of_Schur,rilu->x_C,rilu->b_C);          Paso_Solver_solveRILU(rilu->RILU_of_Schur,rilu->x_C,rilu->b_C);
347          /* b_F=b_F-A_FC*x_C */          /* b_F=b_F-A_FC*x_C */
348          Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(-1.,rilu->A_FC,rilu->x_C,1.,rilu->b_F);          Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,rilu->A_FC,rilu->x_C,1.,rilu->b_F);
349          /* x_F=invA_FF*b_F  */          /* x_F=invA_FF*b_F  */
350          Paso_Solver_applyBlockDiagonalMatrix(n_block,rilu->n_F,rilu->inv_A_FF,rilu->A_FF_pivot,rilu->x_F,rilu->b_F);          Paso_Solver_applyBlockDiagonalMatrix(n_block,rilu->n_F,rilu->inv_A_FF,rilu->A_FF_pivot,rilu->x_F,rilu->b_F);
351          /* x<-[x_F,x_C]     */          /* x<-[x_F,x_C]     */

Legend:
Removed from v.700  
changed lines
  Added in v.1312

  ViewVC Help
Powered by ViewVC 1.1.26