/[escript]/trunk/esys2/finley/src/finleyC/System.c
ViewVC logotype

Diff of /trunk/esys2/finley/src/finleyC/System.c

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

revision 96 by jgs, Tue Oct 26 06:53:54 2004 UTC revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC
# Line 12  Line 12 
12  /**************************************************************/  /**************************************************************/
13    
14  #include "Finley.h"  #include "Finley.h"
 #include "IndexList.h"  
15  #include "System.h"  #include "System.h"
 #include "Mesh.h"  
16    
17  /**************************************************************/  /**************************************************************/
18    
19  /* allocates a SystemMatrix in the CSR/Harwell boeing format from a  /* allocates a SystemMatrix of type type using the given matrix pattern
20     Finley_Mesh. Values are initialized by zero.  */     if type is UNKOWN CSR is used.
21       if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.
22  Finley_SystemMatrix* Finley_SystemMatrix_alloc(Finley_Mesh *mesh,     if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.
23    Finley_SystemMatrixType type, int symmetric,int row_block_size, int reduce_row_order,     Values are initialized by zero.  */
24    int col_block_size, int reduce_col_order) {  
25    double time0,time1;  Finley_SystemMatrix* Finley_SystemMatrix_alloc(Finley_SystemMatrixType type,Finley_SystemMatrixPattern *pattern, int row_block_size, int col_block_size) {
26    maybelong packed_row_block_size, packed_col_block_size, unpacked_row_block_size, unpacked_col_block_size,i,block_size,j,iptr;    double time0;
27    maybelong *rowLabel=NULL,*colLabel=NULL,len_index_list=0;    Finley_SystemMatrix*out=NULL;
28    Finley_SystemMatrix*out;    double *val=NULL;
29    Finley_IndexList* index_list=NULL;    Finley_SystemMatrixType out_type;
30    Finley_ErrorCode=NO_ERROR;    Finley_ErrorCode=NO_ERROR;
31        
32      /* check the matrix type */
33      switch(type) {
34         case CSC:
35            out_type=CSC;
36            break;
37         case CSR:
38            out_type=CSR;
39            break;
40         case CSC_BLK1:
41            out_type=CSC;
42            if (row_block_size!=1 || col_block_size!=1) {
43                Finley_ErrorCode=TYPE_ERROR;
44                sprintf(Finley_ErrorMsg,"convertion of matrix pattern for block size one for logical block size > 1 is not implemented yet");
45                return NULL;
46            }
47            break;
48         case CSR_BLK1:
49            out_type=CSR;
50            if (row_block_size!=1 || col_block_size!=1) {
51                Finley_ErrorCode=TYPE_ERROR;
52                sprintf(Finley_ErrorMsg,"convertion of matrix pattern for block size one for logical block size > 1 is not implemented yet");
53                return NULL;
54            }
55            break;
56         case CSC_SYM:
57         case CSC_BLK1_SYM:
58            out_type=CSC_SYM;
59            Finley_ErrorCode=TYPE_ERROR;
60            sprintf(Finley_ErrorMsg,"convertion of matrix pattern for symmetric CSC is not implemented yet.");
61            return NULL;
62         default:
63            Finley_ErrorCode=TYPE_ERROR;
64            sprintf(Finley_ErrorMsg,"unknown matrix type identifier %d.",type);
65            return NULL;
66      }
67    time0=Finley_timer();    time0=Finley_timer();
   /* the matrix block row_block_size x col_block_size is fully  
   incorporated into the matrix pattern.  this is required for using  
   the SGI scientific library. row_offset is the offset for the row  
   index typically */  
   
   packed_row_block_size=row_block_size;  
   packed_col_block_size=col_block_size;  
   unpacked_row_block_size=1;  
   unpacked_col_block_size=1;  
   
   
68    /*  allocate the return value */    /*  allocate the return value */
69    
70    out=(Finley_SystemMatrix*)MEMALLOC(sizeof(Finley_SystemMatrix));    out=MEMALLOC(1,Finley_SystemMatrix);
71    if (Finley_checkPtr(out)) goto clean;    if (! Finley_checkPtr(out)) {
   
   if (type==UNKNOWN) {  
      out->type=FINLEY_DEFAULT_MATRIX_TYPE;  
   } else {  
      out->type=type;  
   }  
     
   out->total_row_block_size=row_block_size;  
   out->total_col_block_size=col_block_size;  
   
   out->row_block_size=unpacked_row_block_size;  
   out->col_block_size=unpacked_col_block_size;  
   
   block_size=out->row_block_size*out->col_block_size;  
   
 /* ******************************** */  
   if (reduce_col_order) {  
      out->num_cols=packed_col_block_size*mesh->Nodes->reducedNumDegreesOfFreedom;  
      colLabel=mesh->Nodes->reducedDegreeOfFreedom;  
   } else {  
      out->num_cols=packed_col_block_size*mesh->Nodes->numDegreesOfFreedom;  
      colLabel=mesh->Nodes->degreeOfFreedom;  
   }  
       
   if (reduce_row_order) {  
      out->num_rows=packed_row_block_size*mesh->Nodes->reducedNumDegreesOfFreedom;  
      rowLabel=mesh->Nodes->reducedDegreeOfFreedom;  
   } else {  
      out->num_rows=packed_row_block_size*mesh->Nodes->numDegreesOfFreedom;  
      rowLabel=mesh->Nodes->degreeOfFreedom;  
   }  
 /* ******************************** */  
   
   out->symmetric=symmetric;  
   out->index=NULL;  
   out->val=NULL;  
   out->ptr=NULL;  
   out->reference_counter=0;  
   out->lenOfVal=0;  
   out->solve=NULL;  
   out->iterative=NULL;  
   
   /*  initialize the (temporary) list index_list of the colums indices: */  
   switch(out->type) {  
   case CSR:  
     len_index_list=out->num_rows;  
     break;  
   case CSC:  
     len_index_list=out->num_cols;  
     break;  
   default:  
     Finley_ErrorCode=TYPE_ERROR;  
     sprintf(Finley_ErrorMsg,"Unknown matrix type.");  
     goto clean;  
   } /* switch out->type */  
   
   index_list=(Finley_IndexList*) TMPMEMALLOC(len_index_list*sizeof(Finley_IndexList));  
   if (Finley_checkPtr(index_list)) goto clean;  
   
   #pragma omp parallel for private(i) schedule(static)  
   for(i=0;i<len_index_list;i++) {  
        index_list[i].extension=NULL;  
        index_list[i].n=0;  
   }  
   
   /*  insert contributions from element matrices into colums index index_list: */  
   Finley_IndexList_insertElements(index_list,mesh->Elements,  
                                   reduce_row_order,packed_row_block_size,rowLabel,  
                                   reduce_col_order,packed_col_block_size,colLabel,  
                                   symmetric, out->type);  
   Finley_IndexList_insertElements(index_list,mesh->FaceElements,  
                                   reduce_row_order,packed_row_block_size,rowLabel,  
                                   reduce_col_order,packed_col_block_size,colLabel,  
                                   symmetric, out->type);  
   Finley_IndexList_insertElements(index_list,mesh->ContactElements,  
                                   reduce_row_order,packed_row_block_size,rowLabel,  
                                   reduce_col_order,packed_col_block_size,colLabel,  
                                   symmetric, out->type);  
   Finley_IndexList_insertElements(index_list,mesh->Points,  
                                   reduce_row_order,packed_row_block_size,rowLabel,  
                                   reduce_col_order,packed_col_block_size,colLabel,  
                                   symmetric, out->type);  
   if (Finley_ErrorCode!=NO_ERROR) goto clean;  
   
   /* allocate the ptr: */  
   
   out->ptr=(maybelong*)MEMALLOC(((len_index_list)+1)*sizeof(maybelong));  
   if (Finley_checkPtr(out->ptr)) {  
     Finley_SystemMatrix_dealloc(out);  
     goto clean;  
   }  
   #pragma omp parallel for private(i) schedule(static)  
   for(i=0;i<len_index_list;i++) out->ptr[i]=0;  
   out->ptr[len_index_list]=0;  
   
   /* count entries in each column and store to pointer vector: */  
   out->ptr[0]=PTR_OFFSET;  
   /* OMP */  
   for(i=0;i<len_index_list;i++)  
     out->ptr[i+1]=out->ptr[i]+Finley_IndexList_count(&index_list[i]);  
   
   /* allocate index and val: */  
   
   out->index=(maybelong*)MEMALLOC((out->ptr[len_index_list]-PTR_OFFSET)*sizeof(maybelong));  
   out->lenOfVal=(out->ptr[len_index_list]-PTR_OFFSET)*out->row_block_size*out->col_block_size;  
   out->val=(double*)MEMALLOC(out->lenOfVal*sizeof(double));  
   if (Finley_checkPtr(out->index) || Finley_checkPtr(out->val) ) {  
     Finley_SystemMatrix_dealloc(out);  
     goto clean;  
   }  
   #pragma omp parallel firstprivate(out,index_list,len_index_list)  
   {  
72    
73  #pragma omp master       /* is block size 1 enforced ? */
74  {       if (type==CSC_BLK1 || type==CSR_BLK1 || type==CSR_BLK1_SYM || type==CSC_BLK1_SYM) {
75  time1=Finley_timer();          out->row_block_size=1;
76  }          out->col_block_size=1;
77       #pragma omp for private(i,iptr,j) schedule(static)       } else {
78       for (i=0;i< len_index_list;i++) {          out->row_block_size=row_block_size;
79          for (iptr=out->ptr[i]-PTR_OFFSET;iptr<out->ptr[i+1]-PTR_OFFSET; iptr++) {          out->col_block_size=col_block_size;
                out->index[iptr]=0;  
                for (j=0;j<block_size;j++) out->val[iptr*block_size+j]=0.;  
         }  
80       }       }
81  #pragma omp master       if (out_type==CSC || type==CSC_BLK1 ||  type==CSC_SYM || type==CSC_BLK1_SYM ) {
82  {           out->num_cols=pattern->n_index;
83  printf("timing: matrix pattern: instantation: %.4e sec\n",Finley_timer()-time1);           out->num_rows=pattern->n_ptr;
84  }       } else {
85             out->num_rows=pattern->n_index;
86             out->num_cols=pattern->n_ptr;
87         }
88    
89         out->type=out_type;
90         out->logical_row_block_size=row_block_size;
91         out->logical_col_block_size=col_block_size;
92         out->logical_block_size=out->logical_row_block_size*out->logical_block_size;
93         out->block_size=out->row_block_size*out->col_block_size;
94         out->pattern=Finley_SystemMatrixPattern_reference(pattern);
95         out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
96         out->reference_counter=1;
97         out->direct=NULL;  
98         out->iterative=NULL;
99    
      /* change the list of the row indicies into an array:  */  
      /* each index is ordered by increased size  */  
100        
101       #pragma omp for private(i) schedule(static)       /* allocate memory for matrix entries */
102       for(i=0;i<len_index_list;i++) {       val=MEMALLOC(out->len,double);
103         Finley_IndexList_toArray(&index_list[i],&(out->index[out->ptr[i]-PTR_OFFSET]));       if (! Finley_checkPtr(val)) {
104         qsort(&(out->index[out->ptr[i]-PTR_OFFSET]),(int)(out->ptr[i+1]-out->ptr[i]),sizeof(maybelong), Finley_comparIndex);          out->val=val;
105            Finley_SystemMatrix_setValues(out,DBLE(0));
106       }       }
107       }
   }  
   /* out->val is initialized to zero: */  
   
 time1=Finley_timer();  
   Finley_SystemMatrix_setValues(out,DBLE(0));  
 printf("timing: matrix pattern: initialize: %.4e sec\n",Finley_timer()-time1);  
108    /* all done: */    /* all done: */
109    
110      printf("timing: system matrix %.4e sec\n",Finley_timer()-time0);
  clean:  
   if (index_list!=NULL) {  
     #pragma omp parallel for private(i)  
     for(i=0;i<len_index_list;i++) Finley_IndexList_free(index_list[i].extension);  
   }  
   TMPMEMFREE(index_list);  
   
   printf("timing: matrix pattern: %.4e sec\n",Finley_timer()-time0);  
   
111    if (Finley_ErrorCode!=NO_ERROR) {    if (Finley_ErrorCode!=NO_ERROR) {
112        MEMFREE(val);
113        Finley_SystemMatrix_dealloc(out);
114      return NULL;      return NULL;
115    } else {    } else {
116      #ifdef Finley_TRACE      #ifdef Finley_TRACE
117      printf("Finley_SystemMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->num_rows,(long)out->num_cols);      printf("Finley_SystemMatrix_alloc: %ld x %ld system matrix has been allocated.\n",(long)out->num_rows,(long)out->num_cols);
118      #endif      #endif
     out->reference_counter++;  
119      return out;      return out;
120    }    }
121  }  }
122    
123    /* returns a reference to Finley_SystemMatrix in */
124    
125    Finley_SystemMatrix* Finley_SystemMatrix_reference(Finley_SystemMatrix* in) {
126       if (in!=NULL) ++(in->reference_counter);
127       return NULL;
128    }
129    
130  /* deallocates a SystemMatrix: */  /* deallocates a SystemMatrix: */
131    
132  void Finley_SystemMatrix_dealloc(Finley_SystemMatrix* in) {  void Finley_SystemMatrix_dealloc(Finley_SystemMatrix* in) {
133    if (in!=NULL) {    if (in!=NULL) {
134       in->reference_counter--;       in->reference_counter--;
135       if (in->reference_counter<=0) {       if (in->reference_counter<=0) {
         MEMFREE(in->ptr);  
         MEMFREE(in->index);  
136          MEMFREE(in->val);          MEMFREE(in->val);
137          Finley_SystemMatrix_solve_free(in);          Finley_SystemMatrixPattern_dealloc(in->pattern);
138          Finley_SystemMatrix_iterative_free(in);          Finley_SystemMatrix_solve_free(in);
139          MEMFREE(in);          MEMFREE(in);
140          #ifdef Finley_TRACE          #ifdef Finley_TRACE
141          printf("Finley_SystemMatrix_dealloc: system matrix as been deallocated.\n");          printf("Finley_SystemMatrix_dealloc: system matrix as been deallocated.\n");
# Line 233  void Finley_SystemMatrix_dealloc(Finley_ Line 143  void Finley_SystemMatrix_dealloc(Finley_
143       }       }
144     }     }
145  }  }
 /* *************************************************************/  
   
 /*  some routines which help to get the matrix pattern from elements: */  
   
 /*  this routine is used by qsort called in Finley_SystemMatrix_alloc */  
   
 int Finley_comparIndex(const void *index1,const void *index2){  
    maybelong Iindex1,Iindex2;  
    Iindex1=*(maybelong*)index1;  
    Iindex2=*(maybelong*)index2;  
    if (Iindex1<Iindex2) {  
       return -1;  
    } else {  
       if (Iindex1>Iindex2) {  
          return 1;  
       } else {  
          return 0;  
       }  
    }  
 }  
146  /*  /*
147   * $Log$   * $Log$
148   * Revision 1.1  2004/10/26 06:53:57  jgs   * Revision 1.2  2004/12/14 05:39:30  jgs
149   * Initial revision   * *** empty log message ***
150     *
151     * Revision 1.1.1.1.2.3  2004/11/24 01:37:15  gross
152     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
153     *
154     * Revision 1.1.1.1.2.2  2004/11/12 06:58:18  gross
155     * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
156     *
157     * Revision 1.1.1.1.2.1  2004/10/28 22:59:24  gross
158     * finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
159     *
160     * Revision 1.1.1.1  2004/10/26 06:53:57  jgs
161     * initial import of project esys2
162   *   *
163   * Revision 1.1  2004/07/02 04:21:13  gross   * Revision 1.1  2004/07/02 04:21:13  gross
164   * Finley C code has been included   * Finley C code has been included

Legend:
Removed from v.96  
changed lines
  Added in v.97

  ViewVC Help
Powered by ViewVC 1.1.26