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

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

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

revision 2239 by ksteube, Thu Sep 25 23:11:13 2008 UTC revision 2240 by artak, Wed Feb 4 01:17:06 2009 UTC
# Line 25  Line 25 
25  #include "Paso.h"  #include "Paso.h"
26  #include "SparseMatrix.h"  #include "SparseMatrix.h"
27  #include "TRILINOS.h"  #include "TRILINOS.h"
28    #include "mmio.h"
29    
30  /**************************************************************/  /**************************************************************/
31    static void swap( index_t*, index_t*, double*, int, int );
32    static void q_sort( index_t*, index_t*, double*, int, int );
33    /*static void print_entries( index_t*, index_t*, double* );*/
34    
35    static int M, N, nz;
36    
37    
38    /* debug: print the entries */
39    /*
40    void print_entries( index_t *r, index_t *c, double *v )
41    {
42        int i;
43    
44        for( i=0; i<nz; i++ )
45        {
46            printf( "(%ld, %ld) == %e\n", (long)r[i], (long)c[i], v[i] );
47        }
48    }
49    */
50    
51    /* swap function */
52    void swap( index_t *r, index_t *c, double *v, int left, int right )
53    {
54        double v_temp;
55        index_t temp;
56    
57        temp = r[left];
58        r[left] = r[right];
59        r[right] = temp;
60    
61        temp = c[left];
62        c[left] = c[right];
63        c[right] = temp;
64    
65        v_temp = v[left];
66        v[left] = v[right];
67        v[right] = v_temp;
68    }
69    
70    void q_sort( index_t *row, index_t *col, double *val, int begin, int end )
71    {
72        int l, r;
73        index_t pivot, lval;
74    
75        if( end > begin )
76        {
77            pivot = N * row[begin] + col[begin];
78            l = begin + 1;
79            r = end;
80    
81            while( l < r )
82            {
83                lval = N * row[l] + col[l];
84                if( lval < pivot )
85                    l++;
86                else
87                {
88                    r--;
89                    swap( row, col, val, l, r );
90                }
91            }
92            l--;
93            swap( row, col, val, begin, l );
94            q_sort( row, col, val, begin, l );
95            q_sort( row, col, val, r, end );
96        }
97    }
98    
99    
100  /* allocates a SparseMatrix of type type using the given matrix pattern  /* allocates a SparseMatrix of type type using the given matrix pattern
101     if type is UNKOWN CSR is used.     if type is UNKOWN CSR is used.
102     if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.     if CSC or CSC_BLK1 is used pattern has to give the CSC pattern.
103     if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.     if CSR or CSR_BLK1 is used pattern has to give the CSR pattern.
104     Values are initialized by zero.  */     Values are initialized by zero.  */
   
105  Paso_SparseMatrix* Paso_SparseMatrix_alloc(Paso_SparseMatrixType type,Paso_Pattern *pattern, int row_block_size, int col_block_size) {  Paso_SparseMatrix* Paso_SparseMatrix_alloc(Paso_SparseMatrixType type,Paso_Pattern *pattern, int row_block_size, int col_block_size) {
106    
107    double time0;    double time0;
# Line 128  void Paso_SparseMatrix_free(Paso_SparseM Line 196  void Paso_SparseMatrix_free(Paso_SparseM
196       }       }
197     }     }
198  }  }
199    
200    Paso_SparseMatrix* Paso_SparseMatrix_loadMM_toCSR( char *fileName_p )
201    {
202        index_t *col_ind = NULL;
203        index_t *row_ind = NULL;
204        index_t *row_ptr = NULL;
205        double *val = NULL;
206        FILE *fileHandle_p = NULL;
207            Paso_Pattern* mainPattern=NULL;
208        Paso_SparseMatrix *out = NULL;
209        int i, curr_row, scan_ret;
210        MM_typecode matrixCode;
211            Paso_resetError();
212    
213        /* open the file */
214        fileHandle_p = fopen( fileName_p, "r" );
215        if( fileHandle_p == NULL )
216        {
217            Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Cannot read file for reading.");
218            return NULL;
219        }
220    
221        /* process banner */
222        if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
223        {
224            Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Error processing MM banner.");
225            fclose( fileHandle_p );
226            return NULL;
227        }
228        if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
229        {
230    
231            Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
232            fclose( fileHandle_p );
233            return NULL;
234        }
235    
236        /* get matrix size */
237        if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
238        {
239            Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not parse matrix size");
240            fclose( fileHandle_p );
241            return NULL;
242        }
243    
244        /* prepare storage */
245        col_ind = MEMALLOC( nz, index_t );
246        row_ind = MEMALLOC( nz, index_t );
247        val = MEMALLOC( nz, double );
248    
249        row_ptr = MEMALLOC( (M+1), index_t );
250    
251        if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
252        {
253            Paso_setError(MEMORY_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not allocate memory" );
254            fclose( fileHandle_p );
255            return NULL;
256        }
257    
258        /* perform actual read of elements */
259        for( i=0; i<nz; i++ )
260        {
261            scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
262            if (scan_ret!=3)
263            {
264                MEMFREE( val );
265                MEMFREE( row_ind );
266                MEMFREE( col_ind );
267                MEMFREE( row_ptr );
268                fclose(fileHandle_p);
269                return NULL;
270            }
271            row_ind[i]--;
272            col_ind[i]--;
273        }
274        fclose( fileHandle_p );
275    
276        /* sort the entries */
277        q_sort( row_ind, col_ind, val, 0, nz );
278    
279        /* setup row_ptr */
280        curr_row = 0;
281        for( i=0; (i<nz && curr_row<M); curr_row++ )
282        {
283            while( row_ind[i] != curr_row )
284                i++;
285            row_ptr[curr_row] = i;
286        }
287        row_ptr[M] = nz;
288    
289            mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,M,N,row_ptr,col_ind);
290        out  = Paso_SparseMatrix_alloc(MATRIX_FORMAT_DEFAULT, mainPattern, 1, 1);
291        /* copy values and cleanup temps */
292        for( i=0; i<nz; i++ ) out->val[i] = val[i];
293        
294        Paso_Pattern_free(mainPattern);
295        MEMFREE( val );
296        MEMFREE( row_ind );
297        return out;
298    }
299    
300    
301    void Paso_SparseMatrix_saveMM(Paso_SparseMatrix * A_p, char * fileName_p) {
302      FILE * fileHandle_p = NULL;
303      dim_t N,M,i, iptr_ij;
304      MM_typecode matcode;                        
305    
306      if (A_p->col_block_size !=A_p->row_block_size) {
307        Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_saveMM: currently only square blocks are supported.");
308        return;
309      }
310      if (A_p->row_block_size>3) {
311           Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM: currently only block size 3 is supported.\n");
312           return;
313      }
314    
315      if (A_p->type & MATRIX_FORMAT_SYM) {
316        Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support symmetric storage scheme");
317        return;
318      }
319      
320      /* open the file */
321      fileHandle_p = fopen(fileName_p, "w");
322      if (fileHandle_p==NULL) {
323        Paso_setError(IO_ERROR,"file could not be opened for writing");
324        return;
325      }
326    
327      if (A_p->type & MATRIX_FORMAT_CSC) {
328        Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support CSC yet.");
329      } else {
330        mm_initialize_typecode(&matcode);
331        mm_set_matrix(&matcode);
332        mm_set_coordinate(&matcode);
333        mm_set_real(&matcode);
334    
335        N= A_p->numRows;
336        M= A_p->numCols;
337        mm_write_banner(fileHandle_p, matcode);
338        mm_write_mtx_crd_size(fileHandle_p, N, M, A_p->pattern->ptr[N]);
339        if(A_p->row_block_size==1)
340          for (i=0; i<N; i++)
341            for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
342              fprintf(fileHandle_p, "%d %d %25.15e\n", i+1, A_p->pattern->index[iptr_ij]+1, A_p->val[iptr_ij]);
343            }
344    
345        if(A_p->row_block_size==2)
346          for (i=0; i<N; i++) {
347            for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
348              fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1, 2*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*4]);
349              fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1, 2*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*4+2]);
350              fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1+1, 2*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*4+1]);
351              fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1+1, 2*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*4+3]);
352            }
353          }
354    
355        if(A_p->row_block_size==3)
356          for (i=0; i<N; i++) {
357           for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
358              fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1, 3*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9]);
359              fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1, 3*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+3]);
360              fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1, 3*(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+6]);
361              fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1+1, 3*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9+1]);
362              fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1+1, 3*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+4]);
363              fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1+1, 3*(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+7]);
364              fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+2+1, 3*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9+2]);
365              fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+2+1, 3*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+5]);
366              fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+2+1, 3*(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+8]);
367            }
368          }
369        
370      }
371    
372      /* close the file */
373      fclose(fileHandle_p);
374      
375      return;
376    }
377    
Legend:
Removed from v.2239  
changed lines
  Added in v.2240
  ViewVC Help
Powered by ViewVC 1.1.26