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

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

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

revision 1361 by gross, Fri Dec 14 09:26:51 2007 UTC revision 2608 by jfenwick, Tue Aug 18 01:25:18 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 19  Line 18 
18    
19  /**************************************************************/  /**************************************************************/
20    
21  /* Author: gross@access.edu.au */  /* Author: Lutz Gross, l.gross@uq.edu.au */
22    
23  /**************************************************************/  /**************************************************************/
24    
# Line 28  Line 27 
27  /**************************************************************/  /**************************************************************/
28    
29  /* allocates a SystemMatrix of type type using the given matrix pattern  /* allocates a SystemMatrix of type type using the given matrix pattern
30     if type is UNKOWN CSR is used.     Values are initialized by zero.
31     if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.     if patternIsUnrolled and type & MATRIX_FORMAT_BLK1, it is assumed that the pattern is allready unrolled to match the requested block size
32     if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.     and offsets otherwise unrolling and offset adjustment will be performed.
33     Values are initialized by zero.  */  */
34    
35  Paso_SystemMatrix* Paso_SystemMatrix_alloc(Paso_SystemMatrixType type,Paso_SystemMatrixPattern *pattern, int row_block_size, int col_block_size) {  Paso_SystemMatrix* Paso_SystemMatrix_alloc(Paso_SystemMatrixType type,Paso_SystemMatrixPattern *pattern, int row_block_size, int col_block_size,
36      const bool_t patternIsUnrolled) {
37    
   double time0;  
38    Paso_SystemMatrix*out=NULL;    Paso_SystemMatrix*out=NULL;
39    dim_t n_norm,i;    dim_t n_norm,i;
40    Paso_SystemMatrixType pattern_format_out;    Paso_SystemMatrixType pattern_format_out;
41      bool_t unroll=FALSE;
42    
43      pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT;
44    Paso_resetError();    Paso_resetError();
45    time0=Paso_timer();    if (type & MATRIX_FORMAT_SYM) {
46          Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_alloc: Symmetric matrix patterns are not supported.");
47          return NULL;
48      }
49      if (patternIsUnrolled) {
50         if ( ! XNOR(type & MATRIX_FORMAT_OFFSET1, pattern->type & PATTERN_FORMAT_OFFSET1) ) {
51             Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_alloc: requested offset and pattern offset does not match.");
52             return NULL;
53         }
54      }
55      /* do we need to apply unrolling ? */
56      unroll  
57            /* we don't like non-square blocks */
58        =   (row_block_size!=col_block_size)
59            /* or any block size bigger than 3 */
60        ||  (col_block_size>3)
61            /* or if lock size one requested and the block size is not 1 */
62        ||  ((type & MATRIX_FORMAT_BLK1) &&  (col_block_size>1) )
63            /* or the offsets are wrong */
64        ||  ((type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & PATTERN_FORMAT_OFFSET1));
65      
66    out=MEMALLOC(1,Paso_SystemMatrix);    out=MEMALLOC(1,Paso_SystemMatrix);
67    if (! Paso_checkPtr(out)) {      if (! Paso_checkPtr(out)) {  
68       out->type=type;       out->type=type;
# Line 49  Paso_SystemMatrix* Paso_SystemMatrix_all Line 70  Paso_SystemMatrix* Paso_SystemMatrix_all
70       out->row_distribution=NULL;       out->row_distribution=NULL;
71       out->col_distribution=NULL;       out->col_distribution=NULL;
72       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);       out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
73         out->row_coupler=NULL;
74         out->col_coupler=NULL;
75       out->mainBlock=NULL;       out->mainBlock=NULL;
76       out->coupleBlock=NULL;       out->row_coupleBlock=NULL;
77         out->col_coupleBlock=NULL;
78       out->normalizer_is_valid=FALSE;       out->normalizer_is_valid=FALSE;
79       out->normalizer=NULL;       out->normalizer=NULL;
80       out->solver_package=PASO_PASO;         out->solver_package=PASO_PASO;  
81       out->solver=NULL;         out->solver=NULL;  
82       out->trilinos_data=NULL;       out->trilinos_data=NULL;
83       out->reference_counter=1;       out->reference_counter=1;
84         out->logical_row_block_size=row_block_size;
85         out->logical_col_block_size=col_block_size;
86    
87    
      pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT;  
      /* ====== compressed sparse columns === */  
88       if (type & MATRIX_FORMAT_CSC) {       if (type & MATRIX_FORMAT_CSC) {
89          if (type & MATRIX_FORMAT_SYM) {           if (unroll) {
90             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");                 if (patternIsUnrolled) {
91          } else {                    out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
92             if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {                 } else {
93                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);                    out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
94                out->row_block_size=1;                 }
95                out->col_block_size=1;                 out->row_block_size=1;
96             } else {                 out->col_block_size=1;
97             } else {
98                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
99                out->row_block_size=row_block_size;                out->row_block_size=row_block_size;
100                out->col_block_size=col_block_size;                out->col_block_size=col_block_size;
101             }           }
102          }           if (Paso_noError()) {
103          out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);             out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
104          out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);             out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
105             }
106       } else {       } else {
107       /* ====== compressed sparse row === */           if (unroll) {
108          if (type & MATRIX_FORMAT_SYM) {                if (patternIsUnrolled) {
109             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");                    out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
110          } else {                } else {
111             if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3)  {                    out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
112                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);                }
113                out->row_block_size=1;                out->row_block_size=1;
114                out->col_block_size=1;                out->col_block_size=1;
115             } else {           } else {
116                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
117                out->row_block_size=row_block_size;                out->row_block_size=row_block_size;
118                out->col_block_size=col_block_size;                out->col_block_size=col_block_size;
119             }           }
120          }           if (Paso_noError()) {
121          out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);                out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
122          out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);                out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
123             }
124       }       }
125       out->logical_row_block_size=row_block_size;       if (Paso_noError()) {
126       out->logical_col_block_size=col_block_size;          out->block_size=out->row_block_size*out->col_block_size;
127       out->logical_block_size=out->logical_row_block_size*out->logical_block_size;          out->col_coupler=Paso_Coupler_alloc(pattern->col_connector,out->col_block_size);
128       out->block_size=out->row_block_size*out->col_block_size;          out->row_coupler=Paso_Coupler_alloc(pattern->row_connector,out->row_block_size);
129       /* this should be bypassed if trilinos is used */          /* this should be bypassed if trilinos is used */
130       if (type & MATRIX_FORMAT_TRILINOS_CRS) {          if (type & MATRIX_FORMAT_TRILINOS_CRS) {
131          #ifdef TRILINOS             #ifdef TRILINOS
132          out->trilinos_data=Paso_TRILINOS_alloc();             out->trilinos_data=Paso_TRILINOS_alloc();
133          #endif             #endif
      } else {  
         out->solver_package=PASO_PASO;    
         out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size);  
         out->coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->couplePattern,row_block_size,col_block_size);  
         /* allocate memory for matrix entries */  
         if (type & MATRIX_FORMAT_CSC) {  
            n_norm = out->mainBlock->numCols * out->col_block_size;  
134          } else {          } else {
135             n_norm = out->mainBlock->numRows * out->row_block_size;             out->solver_package=PASO_PASO;  
136               out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size,TRUE);
137               out->col_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->col_couplePattern,row_block_size,col_block_size,TRUE);
138               out->row_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->row_couplePattern,row_block_size,col_block_size,TRUE);
139               if (Paso_noError()) {
140                  /* allocate memory for matrix entries */
141                  if (type & MATRIX_FORMAT_CSC) {
142                     n_norm = out->mainBlock->numCols * out->col_block_size;
143                  } else {
144                     n_norm = out->mainBlock->numRows * out->row_block_size;
145                  }
146                  out->normalizer=MEMALLOC(n_norm,double);
147                  out->normalizer_is_valid=FALSE;
148                  if (! Paso_checkPtr(out->normalizer)) {
149                     #pragma omp parallel for private(i) schedule(static)
150                     for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
151                  }
152               }
153          }          }
154          out->normalizer=MEMALLOC(n_norm,double);       }
         out->normalizer_is_valid=FALSE;  
         if (! Paso_checkPtr(out->normalizer)) {  
            #pragma omp parallel for private(i) schedule(static)  
            for (i=0;i<n_norm;++i) out->normalizer[i]=0.;  
        }  
     }  
155    }    }
156    /* all done: */    /* all done: */
157    if (! Paso_noError()) {    if (! Paso_noError()) {
158      Paso_SystemMatrix_free(out);      Paso_SystemMatrix_free(out);
159      return NULL;      return NULL;
160    } else {    } else {
     #ifdef Paso_TRACE  
     printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);  
     printf("Paso_SystemMatrix_alloc: system matrix has been allocated.\n");  
     #endif  
161      return out;      return out;
162    }    }
163  }  }
164    
165  /* returns a reference to Paso_SystemMatrix in */  /* returns a reference to Paso_SystemMatrix in */
166    
167  Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {  Paso_SystemMatrix* Paso_SystemMatrix_getReference(Paso_SystemMatrix* in) {
168     if (in!=NULL) ++(in->reference_counter);     if (in!=NULL) ++(in->reference_counter);
169     return NULL;     return in;
170  }  }
171    
172  /* deallocates a SystemMatrix: */  /* deallocates a SystemMatrix: */
# Line 151  void Paso_SystemMatrix_free(Paso_SystemM Line 179  void Paso_SystemMatrix_free(Paso_SystemM
179          Paso_Distribution_free(in->row_distribution);          Paso_Distribution_free(in->row_distribution);
180          Paso_Distribution_free(in->col_distribution);          Paso_Distribution_free(in->col_distribution);
181          Paso_MPIInfo_free(in->mpi_info);          Paso_MPIInfo_free(in->mpi_info);
182            Paso_Coupler_free(in->row_coupler);
183            Paso_Coupler_free(in->col_coupler);
184          Paso_SparseMatrix_free(in->mainBlock);          Paso_SparseMatrix_free(in->mainBlock);
185          Paso_SparseMatrix_free(in->coupleBlock);          Paso_SparseMatrix_free(in->col_coupleBlock);
186            Paso_SparseMatrix_free(in->row_coupleBlock);
187          MEMFREE(in->normalizer);          MEMFREE(in->normalizer);
188          Paso_solve_free(in);          Paso_solve_free(in);
189          #ifdef TRILINOS          #ifdef TRILINOS
# Line 165  void Paso_SystemMatrix_free(Paso_SystemM Line 196  void Paso_SystemMatrix_free(Paso_SystemM
196       }       }
197     }     }
198  }  }
199  void Paso_SystemMatrix_allocBuffer(Paso_SystemMatrix* A) {  void  Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A,const double* in)
200      Paso_Coupler_allocBuffer(A->pattern->coupler,A->col_block_size);  {
201      Paso_SystemMatrix_startColCollect(A,in);
202    }
203    double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A)
204    {
205     return Paso_SystemMatrix_finishColCollect(A);
206    }
207    
208    void  Paso_SystemMatrix_startColCollect(Paso_SystemMatrix* A,const double* in)
209    {
210      Paso_Coupler_startCollect(A->col_coupler, in);
211  }  }
212  void Paso_SystemMatrix_freeBuffer(Paso_SystemMatrix* A) {  double* Paso_SystemMatrix_finishColCollect(Paso_SystemMatrix* A)
213      Paso_Coupler_freeBuffer(A->pattern->coupler);  {
214     Paso_Coupler_finishCollect(A->col_coupler);
215     return A->col_coupler->recv_buffer;
216  }  }
217  void  Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A, double* in)  void  Paso_SystemMatrix_startRowCollect(Paso_SystemMatrix* A,const double* in)
218  {  {
219    Paso_Coupler_startCollect(A->pattern->coupler, in);    Paso_Coupler_startCollect(A->row_coupler, in);
220  }  }
221  double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A)  double* Paso_SystemMatrix_finishRowCollect(Paso_SystemMatrix* A)
222  {  {
223   Paso_Coupler_finishCollect(A->pattern->coupler);   Paso_Coupler_finishCollect(A->row_coupler);
224   return A->pattern->coupler->recv_buffer;   return A->row_coupler->recv_buffer;
225  }  }
226    
227  dim_t Paso_SystemMatrix_getTotalNumRows(Paso_SystemMatrix* A){  dim_t Paso_SystemMatrix_getTotalNumRows(const Paso_SystemMatrix* A){
228    return A->mainBlock->numRows * A->row_block_size;    return A->mainBlock->numRows * A->row_block_size;
229  }  }
230    
231  dim_t Paso_SystemMatrix_getTotalNumCols(Paso_SystemMatrix* A){  dim_t Paso_SystemMatrix_getTotalNumCols(const Paso_SystemMatrix* A){
232    return A->mainBlock->numCols * A->col_block_size;    return A->mainBlock->numCols * A->col_block_size;
233  }  }
234  dim_t Paso_SystemMatrix_getGlobalNumRows(Paso_SystemMatrix* A) {  dim_t Paso_SystemMatrix_getGlobalNumRows(Paso_SystemMatrix* A) {

Legend:
Removed from v.1361  
changed lines
  Added in v.2608

  ViewVC Help
Powered by ViewVC 1.1.26