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

trunk/esys2/paso/src/SystemMatrix.c revision 150 by jgs, Thu Sep 15 03:44:45 2005 UTC trunk/paso/src/SystemMatrix.c revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1  /* $Id$ */  
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2009 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 6  Line 18 
18    
19  /**************************************************************/  /**************************************************************/
20    
 /* Copyrights by ACcESS Australia 2003, 2004,2005 */  
21  /* Author: gross@access.edu.au */  /* Author: gross@access.edu.au */
22    
23  /**************************************************************/  /**************************************************************/
24    
 #include "Paso.h"  
25  #include "SystemMatrix.h"  #include "SystemMatrix.h"
26    
27  /**************************************************************/  /**************************************************************/
# Line 23  Line 33 
33     Values are initialized by zero.  */     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    
37    double time0;    double time0;
38    Paso_SystemMatrix*out=NULL;    Paso_SystemMatrix*out=NULL;
39      dim_t n_norm,i;
40      Paso_SystemMatrixType pattern_format_out;
41    
42    Paso_resetError();    Paso_resetError();
43    time0=Paso_timer();    time0=Paso_timer();
   dim_t n_norm,i;  
44    out=MEMALLOC(1,Paso_SystemMatrix);    out=MEMALLOC(1,Paso_SystemMatrix);
45    if (! Paso_checkPtr(out)) {      if (! Paso_checkPtr(out)) {  
46         out->type=type;
47       out->pattern=NULL;         out->pattern=NULL;  
48       out->direct=NULL;         out->row_distribution=NULL;
49       out->iterative=NULL;       out->col_distribution=NULL;
50       out->val=NULL;         out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
51         out->row_coupler=NULL;
52         out->col_coupler=NULL;
53         out->mainBlock=NULL;
54         out->row_coupleBlock=NULL;
55         out->col_coupleBlock=NULL;
56         out->normalizer_is_valid=FALSE;
57         out->normalizer=NULL;
58         out->solver_package=PASO_PASO;  
59         out->solver=NULL;  
60         out->trilinos_data=NULL;
61       out->reference_counter=1;       out->reference_counter=1;
62       /* check the matrix type */  
63       switch(type) {  
64          case CSC:       pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT;
65            out->type=CSC;       /* ====== compressed sparse columns === */
66            if (row_block_size!=col_block_size || col_block_size>3) {       if (type & MATRIX_FORMAT_CSC) {
67               out->row_block_size=1;          if (type & MATRIX_FORMAT_SYM) {
68               out->col_block_size=1;             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");
69               out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,col_block_size,row_block_size);          } else {
70            } else {             if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
71               out->pattern=Paso_SystemMatrixPattern_reference(pattern);                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
              out->row_block_size=row_block_size;  
              out->col_block_size=col_block_size;  
           }  
           break;  
         case CSR:  
            out->type=CSR;  
            if (row_block_size!=col_block_size || col_block_size>3) {  
72                out->row_block_size=1;                out->row_block_size=1;
73                out->col_block_size=1;                out->col_block_size=1;
74                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,row_block_size,col_block_size);             } else {
75             } else {                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
               out->pattern=Paso_SystemMatrixPattern_reference(pattern);  
76                out->row_block_size=row_block_size;                out->row_block_size=row_block_size;
77                out->col_block_size=col_block_size;                out->col_block_size=col_block_size;
78            }             }
79            break;          }
80          case CSC_BLK1:          out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
81            out->type=CSC;          out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
           out->row_block_size=1;  
           out->col_block_size=1;  
           if (row_block_size==1 && col_block_size==1) {  
               out->pattern=Paso_SystemMatrixPattern_reference(pattern);  
           } else {  
              out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,col_block_size,row_block_size);  
           }  
           break;  
         case CSR_BLK1:  
           out->type=CSR;  
           out->row_block_size=1;  
           out->col_block_size=1;  
           if (row_block_size==1 && col_block_size==1) {  
               out->pattern=Paso_SystemMatrixPattern_reference(pattern);  
           } else {  
              out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,row_block_size,col_block_size);  
           }  
           break;  
         case CSC_SYM:  
         case CSC_BLK1_SYM:  
           out->type=CSC_SYM;  
           Paso_setError(TYPE_ERROR,"convertion of matrix pattern for symmetric CSC is not implemented yet.");  
           return NULL;  
         default:  
           Paso_setError(TYPE_ERROR,"unknown matrix type identifier.");  
           return NULL;  
      }  
      if (out->type==CSC || out->type==CSC_SYM ) {  
          out->num_rows=out->pattern->n_index;  
          out->num_cols=out->pattern->n_ptr;  
          n_norm = out->num_cols * out->col_block_size;  
   
82       } else {       } else {
83           out->num_rows=out->pattern->n_ptr;       /* ====== compressed sparse row === */
84           out->num_cols=out->pattern->n_index;          if (type & MATRIX_FORMAT_SYM) {
85           n_norm = out->num_rows * out->row_block_size;             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");
86       }          } else {
87               if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3)  {
88                  out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
89                  out->row_block_size=1;
90                  out->col_block_size=1;
91               } else {
92                  out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
93                  out->row_block_size=row_block_size;
94                  out->col_block_size=col_block_size;
95               }
96            }
97            out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
98            out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
99         }
100       out->logical_row_block_size=row_block_size;       out->logical_row_block_size=row_block_size;
101       out->logical_col_block_size=col_block_size;       out->logical_col_block_size=col_block_size;
102       out->logical_block_size=out->logical_row_block_size*out->logical_block_size;       out->logical_block_size=out->logical_row_block_size*out->logical_col_block_size;
103       out->block_size=out->row_block_size*out->col_block_size;       out->block_size=out->row_block_size*out->col_block_size;
104       out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);       out->col_coupler=Paso_Coupler_alloc(pattern->col_connector,out->col_block_size);
105       /* allocate memory for matrix entries */       out->row_coupler=Paso_Coupler_alloc(pattern->row_connector,out->row_block_size);
106       out->val=MEMALLOC(out->len,double);       /* this should be bypassed if trilinos is used */
107       out->normalizer=MEMALLOC(n_norm,double);       if (type & MATRIX_FORMAT_TRILINOS_CRS) {
108       out->normalizer_is_valid=FALSE;          #ifdef TRILINOS
109       if (! Paso_checkPtr(out->val)) {          out->trilinos_data=Paso_TRILINOS_alloc();
110          Paso_SystemMatrix_setValues(out,DBLE(0));          #endif
111       }       } else {
112       if (! Paso_checkPtr(out->normalizer)) {          out->solver_package=PASO_PASO;  
113           #pragma omp parallel for private(i) schedule(static)          out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size);
114           for (i=0;i<n_norm;++i) out->normalizer[i]=0.;          out->col_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->col_couplePattern,row_block_size,col_block_size);
115       }          out->row_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->row_couplePattern,row_block_size,col_block_size);
116            /* allocate memory for matrix entries */
117            if (type & MATRIX_FORMAT_CSC) {
118               n_norm = out->mainBlock->numCols * out->col_block_size;
119            } else {
120               n_norm = out->mainBlock->numRows * out->row_block_size;
121            }
122            out->normalizer=MEMALLOC(n_norm,double);
123            out->normalizer_is_valid=FALSE;
124            if (! Paso_checkPtr(out->normalizer)) {
125               #pragma omp parallel for private(i) schedule(static)
126               for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
127           }
128        }
129    }    }
130    /* all done: */    /* all done: */
131    if (! Paso_noError()) {    if (! Paso_noError()) {
132      Paso_SystemMatrix_dealloc(out);      Paso_SystemMatrix_free(out);
133      return NULL;      return NULL;
134    } else {    } else {
135      #ifdef Paso_TRACE      #ifdef Paso_TRACE
136      printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);      printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
137      printf("Paso_SystemMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->num_rows,(long)out->num_cols);      printf("Paso_SystemMatrix_alloc: system matrix has been allocated.\n");
138      #endif      #endif
139      return out;      return out;
140    }    }
# Line 134  Paso_SystemMatrix* Paso_SystemMatrix_all Line 144  Paso_SystemMatrix* Paso_SystemMatrix_all
144    
145  Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {  Paso_SystemMatrix* Paso_SystemMatrix_reference(Paso_SystemMatrix* in) {
146     if (in!=NULL) ++(in->reference_counter);     if (in!=NULL) ++(in->reference_counter);
147     return NULL;     return in;
148  }  }
149    
150  /* deallocates a SystemMatrix: */  /* deallocates a SystemMatrix: */
151    
152  void Paso_SystemMatrix_dealloc(Paso_SystemMatrix* in) {  void Paso_SystemMatrix_free(Paso_SystemMatrix* in) {
153    if (in!=NULL) {    if (in!=NULL) {
154       in->reference_counter--;       in->reference_counter--;
155       if (in->reference_counter<=0) {       if (in->reference_counter<=0) {
156          MEMFREE(in->val);          Paso_SystemMatrixPattern_free(in->pattern);
157            Paso_Distribution_free(in->row_distribution);
158            Paso_Distribution_free(in->col_distribution);
159            Paso_MPIInfo_free(in->mpi_info);
160            Paso_Coupler_free(in->row_coupler);
161            Paso_Coupler_free(in->col_coupler);
162            Paso_SparseMatrix_free(in->mainBlock);
163            Paso_SparseMatrix_free(in->col_coupleBlock);
164            Paso_SparseMatrix_free(in->row_coupleBlock);
165          MEMFREE(in->normalizer);          MEMFREE(in->normalizer);
         Paso_SystemMatrixPattern_dealloc(in->pattern);  
166          Paso_solve_free(in);          Paso_solve_free(in);
167            #ifdef TRILINOS
168            Paso_TRILINOS_free(in->trilinos_data);
169            #endif
170          MEMFREE(in);          MEMFREE(in);
171          #ifdef Paso_TRACE          #ifdef Paso_TRACE
172          printf("Paso_SystemMatrix_dealloc: system matrix as been deallocated.\n");          printf("Paso_SystemMatrix_free: system matrix as been deallocated.\n");
173          #endif          #endif
174       }       }
175     }     }
176  }  }
177  /*  void  Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A,const double* in)
178   * $Log$  {
179   * Revision 1.2  2005/09/15 03:44:38  jgs    Paso_SystemMatrix_startColCollect(A,in);
180   * Merge of development branch dev-02 back to main trunk on 2005-09-15  }
181   *  double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A)
182   * Revision 1.1.2.2  2005/09/07 00:59:08  gross  {
183   * some inconsistent renaming fixed to make the linking work.   return Paso_SystemMatrix_finishColCollect(A);
184   *  }
185   * Revision 1.1.2.1  2005/09/05 06:29:47  gross  
186   * These files have been extracted from finley to define a stand alone libray for iterative  void  Paso_SystemMatrix_startColCollect(Paso_SystemMatrix* A,const double* in)
187   * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but  {
188   * has not been tested yet.    Paso_Coupler_startCollect(A->col_coupler, in);
189   *  }
190   *  double* Paso_SystemMatrix_finishColCollect(Paso_SystemMatrix* A)
191   */  {
192     Paso_Coupler_finishCollect(A->col_coupler);
193     return A->col_coupler->recv_buffer;
194    }
195    void  Paso_SystemMatrix_startRowCollect(Paso_SystemMatrix* A,const double* in)
196    {
197      Paso_Coupler_startCollect(A->row_coupler, in);
198    }
199    double* Paso_SystemMatrix_finishRowCollect(Paso_SystemMatrix* A)
200    {
201     Paso_Coupler_finishCollect(A->row_coupler);
202     return A->row_coupler->recv_buffer;
203    }
204    
205    dim_t Paso_SystemMatrix_getTotalNumRows(const Paso_SystemMatrix* A){
206      return A->mainBlock->numRows * A->row_block_size;
207    }
208    
209    dim_t Paso_SystemMatrix_getTotalNumCols(const Paso_SystemMatrix* A){
210      return A->mainBlock->numCols * A->col_block_size;
211    }
212    dim_t Paso_SystemMatrix_getGlobalNumRows(Paso_SystemMatrix* A) {
213      if (A->type & MATRIX_FORMAT_CSC) {
214          return  Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
215      }  else {
216          return  Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
217      }
218    }
219    dim_t Paso_SystemMatrix_getGlobalNumCols(Paso_SystemMatrix* A) {
220      if (A->type & MATRIX_FORMAT_CSC) {
221          return  Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
222      }  else {
223          return  Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
224      }
225    
226    }
227    dim_t Paso_SystemMatrix_getNumOutput(Paso_SystemMatrix* A) {
228       return Paso_SystemMatrixPattern_getNumOutput(A->pattern);
229    }
230    

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

  ViewVC Help
Powered by ViewVC 1.1.26