/[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 2550 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 2551 by gross, Thu Jul 23 09:19:15 2009 UTC
# Line 27  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 ( (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      
64    out=MEMALLOC(1,Paso_SystemMatrix);    out=MEMALLOC(1,Paso_SystemMatrix);
65    if (! Paso_checkPtr(out)) {      if (! Paso_checkPtr(out)) {  
66       out->type=type;       out->type=type;
# Line 59  Paso_SystemMatrix* Paso_SystemMatrix_all Line 79  Paso_SystemMatrix* Paso_SystemMatrix_all
79       out->solver=NULL;         out->solver=NULL;  
80       out->trilinos_data=NULL;       out->trilinos_data=NULL;
81       out->reference_counter=1;       out->reference_counter=1;
82         out->logical_row_block_size=row_block_size;
83         out->logical_col_block_size=col_block_size;
84    
85    
      pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT;  
      /* ====== compressed sparse columns === */  
86       if (type & MATRIX_FORMAT_CSC) {       if (type & MATRIX_FORMAT_CSC) {
87          if (type & MATRIX_FORMAT_SYM) {           if (unroll) {
88             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");                 if (patternIsUnrolled) {
89          } else {                    out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
90             if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {                 } else {
91                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);
92                out->row_block_size=1;                 }
93                out->col_block_size=1;                 out->row_block_size=1;
94             } else {                 out->col_block_size=1;
95             } else {
96                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
97                out->row_block_size=row_block_size;                out->row_block_size=row_block_size;
98                out->col_block_size=col_block_size;                out->col_block_size=col_block_size;
99             }           }
100          }           if (Paso_noError()) {
101          out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);             out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
102          out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);             out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
103             }
104       } else {       } else {
105       /* ====== compressed sparse row === */           if (unroll) {
106          if (type & MATRIX_FORMAT_SYM) {                if (patternIsUnrolled) {
107             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");                    out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
108          } else {                } else {
109             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);
110                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);                }
111                out->row_block_size=1;                out->row_block_size=1;
112                out->col_block_size=1;                out->col_block_size=1;
113             } else {           } else {
114                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
115                out->row_block_size=row_block_size;                out->row_block_size=row_block_size;
116                out->col_block_size=col_block_size;                out->col_block_size=col_block_size;
117             }           }
118          }           if (Paso_noError()) {
119          out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);                out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
120          out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);                out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
121             }
122       }       }
123       out->logical_row_block_size=row_block_size;       if (Paso_noError()) {
124       out->logical_col_block_size=col_block_size;          out->block_size=out->row_block_size*out->col_block_size;
125       out->logical_block_size=out->logical_row_block_size*out->logical_col_block_size;          out->col_coupler=Paso_Coupler_alloc(pattern->col_connector,out->col_block_size);
126       out->block_size=out->row_block_size*out->col_block_size;          out->row_coupler=Paso_Coupler_alloc(pattern->row_connector,out->row_block_size);
127       out->col_coupler=Paso_Coupler_alloc(pattern->col_connector,out->col_block_size);          /* this should be bypassed if trilinos is used */
128       out->row_coupler=Paso_Coupler_alloc(pattern->row_connector,out->row_block_size);          if (type & MATRIX_FORMAT_TRILINOS_CRS) {
129       /* this should be bypassed if trilinos is used */             #ifdef TRILINOS
130       if (type & MATRIX_FORMAT_TRILINOS_CRS) {             out->trilinos_data=Paso_TRILINOS_alloc();
131          #ifdef TRILINOS             #endif
         out->trilinos_data=Paso_TRILINOS_alloc();  
         #endif  
      } else {  
         out->solver_package=PASO_PASO;    
         out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size);  
         out->col_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->col_couplePattern,row_block_size,col_block_size);  
         out->row_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->row_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;  
132          } else {          } else {
133             n_norm = out->mainBlock->numRows * out->row_block_size;             out->solver_package=PASO_PASO;  
134               out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size,TRUE);
135               out->col_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->col_couplePattern,row_block_size,col_block_size,TRUE);
136               out->row_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->row_couplePattern,row_block_size,col_block_size,TRUE);
137               if (Paso_noError()) {
138                  /* allocate memory for matrix entries */
139                  if (type & MATRIX_FORMAT_CSC) {
140                     n_norm = out->mainBlock->numCols * out->col_block_size;
141                  } else {
142                     n_norm = out->mainBlock->numRows * out->row_block_size;
143                  }
144                  out->normalizer=MEMALLOC(n_norm,double);
145                  out->normalizer_is_valid=FALSE;
146                  if (! Paso_checkPtr(out->normalizer)) {
147                     #pragma omp parallel for private(i) schedule(static)
148                     for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
149                  }
150               }
151          }          }
152          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.;  
        }  
     }  
153    }    }
154    /* all done: */    /* all done: */
155    if (! Paso_noError()) {    if (! Paso_noError()) {
156      Paso_SystemMatrix_free(out);      Paso_SystemMatrix_free(out);
157      return NULL;      return NULL;
158    } 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  
159      return out;      return out;
160    }    }
161  }  }
162    
163  /* returns a reference to Paso_SystemMatrix in */  /* returns a reference to Paso_SystemMatrix in */
164    
165  Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {  Paso_SystemMatrix* Paso_SystemMatrix_getReference(Paso_SystemMatrix* in) {
166     if (in!=NULL) ++(in->reference_counter);     if (in!=NULL) ++(in->reference_counter);
167     return in;     return in;
168  }  }

Legend:
Removed from v.2550  
changed lines
  Added in v.2551

  ViewVC Help
Powered by ViewVC 1.1.26