/[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 1118 by gross, Tue Apr 24 08:55:04 2007 UTC revision 1361 by gross, Fri Dec 14 09:26:51 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    
 /* Copyrights by ACcESS Australia 2003, 2004,2005 */  
22  /* Author: gross@access.edu.au */  /* Author: gross@access.edu.au */
23    
24  /**************************************************************/  /**************************************************************/
25    
 #include "Paso.h"  
26  #include "SystemMatrix.h"  #include "SystemMatrix.h"
27    
28  /**************************************************************/  /**************************************************************/
# Line 39  Paso_SystemMatrix* Paso_SystemMatrix_all Line 38  Paso_SystemMatrix* Paso_SystemMatrix_all
38    double time0;    double time0;
39    Paso_SystemMatrix*out=NULL;    Paso_SystemMatrix*out=NULL;
40    dim_t n_norm,i;    dim_t n_norm,i;
41      Paso_SystemMatrixType pattern_format_out;
42    
43    Paso_resetError();    Paso_resetError();
44    time0=Paso_timer();    time0=Paso_timer();
45    out=MEMALLOC(1,Paso_SystemMatrix);    out=MEMALLOC(1,Paso_SystemMatrix);
46    if (! Paso_checkPtr(out)) {      if (! Paso_checkPtr(out)) {  
47         out->type=type;
48       out->pattern=NULL;         out->pattern=NULL;  
49         out->row_distribution=NULL;
50         out->col_distribution=NULL;
51         out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
52         out->mainBlock=NULL;
53         out->coupleBlock=NULL;
54         out->normalizer_is_valid=FALSE;
55         out->normalizer=NULL;
56       out->solver_package=PASO_PASO;         out->solver_package=PASO_PASO;  
57       out->solver=NULL;         out->solver=NULL;  
58       out->val=NULL;         out->trilinos_data=NULL;
59       out->reference_counter=1;       out->reference_counter=1;
60       out->type=type;  
61         pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT;
62         /* ====== compressed sparse columns === */
63       if (type & MATRIX_FORMAT_CSC) {       if (type & MATRIX_FORMAT_CSC) {
64          if (type & MATRIX_FORMAT_SYM) {          if (type & MATRIX_FORMAT_SYM) {
65             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSC is not implemented yet.");
            return NULL;  
66          } else {          } else {
67             if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {             if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3) {
68                if (type & MATRIX_FORMAT_OFFSET1) {                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
                   out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,col_block_size,row_block_size);  
               } else {  
                   out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,col_block_size,row_block_size);  
               }  
69                out->row_block_size=1;                out->row_block_size=1;
70                out->col_block_size=1;                out->col_block_size=1;
71             } else {             } else {
72                if ( (type & MATRIX_FORMAT_OFFSET1) ==(pattern->type & PATTERN_FORMAT_OFFSET1)) {                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
                   out->pattern=Paso_SystemMatrixPattern_reference(pattern);  
               } else {  
                   out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT,1,1);  
               }  
73                out->row_block_size=row_block_size;                out->row_block_size=row_block_size;
74                out->col_block_size=col_block_size;                out->col_block_size=col_block_size;
75             }             }
76          }          }
77          out->num_rows=out->pattern->n_index;          out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
78          out->num_cols=out->pattern->n_ptr;          out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
         n_norm = out->num_cols * out->col_block_size;  
79       } else {       } else {
80         /* ====== compressed sparse row === */
81          if (type & MATRIX_FORMAT_SYM) {          if (type & MATRIX_FORMAT_SYM) {
82             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");             Paso_setError(TYPE_ERROR,"Generation of matrix pattern for symmetric CSR is not implemented yet.");
            return NULL;  
83          } else {          } else {
84             if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3)  {             if ((type & MATRIX_FORMAT_BLK1) || row_block_size!=col_block_size || col_block_size>3)  {
85                if (type & MATRIX_FORMAT_OFFSET1) {                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
                   out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_OFFSET1,row_block_size,col_block_size);  
               } else {  
                   out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,PATTERN_FORMAT_DEFAULT,row_block_size,col_block_size);  
               }  
86                out->row_block_size=1;                out->row_block_size=1;
87                out->col_block_size=1;                out->col_block_size=1;
88             } else {             } else {
89                if ((type & MATRIX_FORMAT_OFFSET1)==(pattern->type & PATTERN_FORMAT_OFFSET1)) {                out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
                   out->pattern=Paso_SystemMatrixPattern_reference(pattern);  
               } else {  
                   out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,(type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT,1,1);  
               }  
90                out->row_block_size=row_block_size;                out->row_block_size=row_block_size;
91                out->col_block_size=col_block_size;                out->col_block_size=col_block_size;
92             }             }
93          }          }
94          out->num_rows=out->pattern->n_index;          out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
95          out->num_cols=out->pattern->n_ptr;          out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
         n_norm = out->num_cols * out->col_block_size;  
         out->num_rows=out->pattern->n_ptr;  
         out->num_cols=out->pattern->n_index;  
         n_norm = out->num_rows * out->row_block_size;  
96       }       }
97       out->logical_row_block_size=row_block_size;       out->logical_row_block_size=row_block_size;
98       out->logical_col_block_size=col_block_size;       out->logical_col_block_size=col_block_size;
99       out->logical_block_size=out->logical_row_block_size*out->logical_block_size;       out->logical_block_size=out->logical_row_block_size*out->logical_block_size;
100       out->block_size=out->row_block_size*out->col_block_size;       out->block_size=out->row_block_size*out->col_block_size;
101       out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);       /* this should be bypassed if trilinos is used */
102       /* allocate memory for matrix entries */       if (type & MATRIX_FORMAT_TRILINOS_CRS) {
103       out->val=MEMALLOC(out->len,double);          #ifdef TRILINOS
104       out->normalizer=MEMALLOC(n_norm,double);          out->trilinos_data=Paso_TRILINOS_alloc();
105       out->normalizer_is_valid=FALSE;          #endif
106       if (! Paso_checkPtr(out->val)) {       } else {
107          Paso_SystemMatrix_setValues(out,DBLE(0));          out->solver_package=PASO_PASO;  
108       }          out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size);
109       if (! Paso_checkPtr(out->normalizer)) {          out->coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->couplePattern,row_block_size,col_block_size);
110           #pragma omp parallel for private(i) schedule(static)          /* allocate memory for matrix entries */
111           for (i=0;i<n_norm;++i) out->normalizer[i]=0.;          if (type & MATRIX_FORMAT_CSC) {
112       }             n_norm = out->mainBlock->numCols * out->col_block_size;
113            } else {
114               n_norm = out->mainBlock->numRows * out->row_block_size;
115            }
116            out->normalizer=MEMALLOC(n_norm,double);
117            out->normalizer_is_valid=FALSE;
118            if (! Paso_checkPtr(out->normalizer)) {
119               #pragma omp parallel for private(i) schedule(static)
120               for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
121           }
122        }
123    }    }
124    /* all done: */    /* all done: */
125    if (! Paso_noError()) {    if (! Paso_noError()) {
126      Paso_SystemMatrix_dealloc(out);      Paso_SystemMatrix_free(out);
127      return NULL;      return NULL;
128    } else {    } else {
129      #ifdef Paso_TRACE      #ifdef Paso_TRACE
130      printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);      printf("timing: system matrix %.4e sec\n",Paso_timer()-time0);
131      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");
132      #endif      #endif
133      return out;      return out;
134    }    }
# Line 145  Paso_SystemMatrix* Paso_SystemMatrix_ref Line 143  Paso_SystemMatrix* Paso_SystemMatrix_ref
143    
144  /* deallocates a SystemMatrix: */  /* deallocates a SystemMatrix: */
145    
146  void Paso_SystemMatrix_dealloc(Paso_SystemMatrix* in) {  void Paso_SystemMatrix_free(Paso_SystemMatrix* in) {
147    if (in!=NULL) {    if (in!=NULL) {
      printf("Paso_SystemMatrix_dealloc %d: %d\n",in,in->reference_counter);  
148       in->reference_counter--;       in->reference_counter--;
149       if (in->reference_counter<=0) {       if (in->reference_counter<=0) {
150          printf("Paso_SystemMatrix_dealloc %d: system matrix as been deallocated.\n",in);          Paso_SystemMatrixPattern_free(in->pattern);
151          MEMFREE(in->val);          Paso_Distribution_free(in->row_distribution);
152            Paso_Distribution_free(in->col_distribution);
153            Paso_MPIInfo_free(in->mpi_info);
154            Paso_SparseMatrix_free(in->mainBlock);
155            Paso_SparseMatrix_free(in->coupleBlock);
156          MEMFREE(in->normalizer);          MEMFREE(in->normalizer);
         Paso_SystemMatrixPattern_dealloc(in->pattern);  
157          Paso_solve_free(in);          Paso_solve_free(in);
158            #ifdef TRILINOS
159            Paso_TRILINOS_free(in->trilinos_data);
160            #endif
161          MEMFREE(in);          MEMFREE(in);
162          #ifdef Paso_TRACE          #ifdef Paso_TRACE
163          printf("Paso_SystemMatrix_dealloc: system matrix as been deallocated.\n");          printf("Paso_SystemMatrix_free: system matrix as been deallocated.\n");
164          #endif          #endif
165       }       }
166     }     }
167  }  }
168  /*  void Paso_SystemMatrix_allocBuffer(Paso_SystemMatrix* A) {
169   * $Log$      Paso_Coupler_allocBuffer(A->pattern->coupler,A->col_block_size);
170   * Revision 1.2  2005/09/15 03:44:38  jgs  }
171   * Merge of development branch dev-02 back to main trunk on 2005-09-15  void Paso_SystemMatrix_freeBuffer(Paso_SystemMatrix* A) {
172   *      Paso_Coupler_freeBuffer(A->pattern->coupler);
173   * Revision 1.1.2.2  2005/09/07 00:59:08  gross  }
174   * some inconsistent renaming fixed to make the linking work.  void  Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A, double* in)
175   *  {
176   * Revision 1.1.2.1  2005/09/05 06:29:47  gross    Paso_Coupler_startCollect(A->pattern->coupler, in);
177   * These files have been extracted from finley to define a stand alone libray for iterative  }
178   * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but  double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A)
179   * has not been tested yet.  {
180   *   Paso_Coupler_finishCollect(A->pattern->coupler);
181   *   return A->pattern->coupler->recv_buffer;
182   */  }
183    
184    dim_t Paso_SystemMatrix_getTotalNumRows(Paso_SystemMatrix* A){
185      return A->mainBlock->numRows * A->row_block_size;
186    }
187    
188    dim_t Paso_SystemMatrix_getTotalNumCols(Paso_SystemMatrix* A){
189      return A->mainBlock->numCols * A->col_block_size;
190    }
191    dim_t Paso_SystemMatrix_getGlobalNumRows(Paso_SystemMatrix* A) {
192      if (A->type & MATRIX_FORMAT_CSC) {
193          return  Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
194      }  else {
195          return  Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
196      }
197    }
198    dim_t Paso_SystemMatrix_getGlobalNumCols(Paso_SystemMatrix* A) {
199      if (A->type & MATRIX_FORMAT_CSC) {
200          return  Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
201      }  else {
202          return  Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
203      }
204    
205    }
206    dim_t Paso_SystemMatrix_getNumOutput(Paso_SystemMatrix* A) {
207       return Paso_SystemMatrixPattern_getNumOutput(A->pattern);
208    }
209    

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

  ViewVC Help
Powered by ViewVC 1.1.26