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

revision 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1
1
/* \$Id: SparseMatrix.c 1306 2007-09-18 05:51:09Z ksteube \$ */

2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
11   *  *
12   *******************************************************/  *******************************************************/
13
14
15  /**************************************************************/  /**************************************************************/
16
# Line 26  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;
108    Paso_SparseMatrix*out=NULL;    Paso_SparseMatrix*out=NULL;
dim_t n_norm,i;
109    Paso_SparseMatrixType pattern_format_out;    Paso_SparseMatrixType pattern_format_out;
110
111    Paso_resetError();    Paso_resetError();
# Line 50  Paso_SparseMatrix* Paso_SparseMatrix_all Line 116  Paso_SparseMatrix* Paso_SparseMatrix_all
116       out->val=NULL;         out->val=NULL;
117       out->reference_counter=1;       out->reference_counter=1;
118       out->type=type;       out->type=type;
119         out->solver=NULL;
120
121       pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT;       pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1:  PATTERN_FORMAT_DEFAULT;
122       /* ====== compressed sparse columns === */       /* ====== compressed sparse columns === */
# Line 130  void Paso_SparseMatrix_free(Paso_SparseM Line 197  void Paso_SparseMatrix_free(Paso_SparseM
197       }       }
198     }     }
199  }  }
200
201    Paso_SparseMatrix* Paso_SparseMatrix_loadMM_toCSR( char *fileName_p )
202    {
203        index_t *col_ind = NULL;
204        index_t *row_ind = NULL;
205        index_t *row_ptr = NULL;
206        double *val = NULL;
207        FILE *fileHandle_p = NULL;
208            Paso_Pattern* mainPattern=NULL;
209        Paso_SparseMatrix *out = NULL;
210        int i, curr_row, scan_ret;
211        MM_typecode matrixCode;
212            Paso_resetError();
213
214        /* open the file */
215        fileHandle_p = fopen( fileName_p, "r" );
216        if( fileHandle_p == NULL )
217        {
219            return NULL;
220        }
221
222        /* process banner */
223        if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
224        {
225            Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Error processing MM banner.");
226            fclose( fileHandle_p );
227            return NULL;
228        }
229        if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
230        {
231
232            Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
233            fclose( fileHandle_p );
234            return NULL;
235        }
236
237        /* get matrix size */
238        if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
239        {
240            Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not parse matrix size");
241            fclose( fileHandle_p );
242            return NULL;
243        }
244
245        /* prepare storage */
246        col_ind = MEMALLOC( nz, index_t );
247        row_ind = MEMALLOC( nz, index_t );
248        val = MEMALLOC( nz, double );
249
250        row_ptr = MEMALLOC( (M+1), index_t );
251
252        if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
253        {
254            Paso_setError(MEMORY_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not allocate memory" );
255            fclose( fileHandle_p );
256            return NULL;
257        }
258
259        /* perform actual read of elements */
260        for( i=0; i<nz; i++ )
261        {
262            scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
263            if (scan_ret!=3)
264            {
265                MEMFREE( val );
266                MEMFREE( row_ind );
267                MEMFREE( col_ind );
268                MEMFREE( row_ptr );
269                fclose(fileHandle_p);
270                return NULL;
271            }
272            row_ind[i]--;
273            col_ind[i]--;
274        }
275        fclose( fileHandle_p );
276
277        /* sort the entries */
278        q_sort( row_ind, col_ind, val, 0, nz );
279
280        /* setup row_ptr */
281        curr_row = 0;
282        for( i=0; (i<nz && curr_row<M); curr_row++ )
283        {
284            while( row_ind[i] != curr_row )
285                i++;
286            row_ptr[curr_row] = i;
287        }
288        row_ptr[M] = nz;
289
290            mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,M,N,row_ptr,col_ind);
291        out  = Paso_SparseMatrix_alloc(MATRIX_FORMAT_DEFAULT, mainPattern, 1, 1);
292        /* copy values and cleanup temps */
293        for( i=0; i<nz; i++ ) out->val[i] = val[i];
294
295        Paso_Pattern_free(mainPattern);
296        MEMFREE( val );
297        MEMFREE( row_ind );
298        return out;
299    }
300
301
302    void Paso_SparseMatrix_saveMM(Paso_SparseMatrix * A_p, char * fileName_p) {
303      FILE * fileHandle_p = NULL;
304      dim_t N,M,i, iptr_ij;
305      MM_typecode matcode;
306
307      if (A_p->col_block_size !=A_p->row_block_size) {
308        Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_saveMM: currently only square blocks are supported.");
309        return;
310      }
311      if (A_p->row_block_size>3) {
312           Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM: currently only block size 3 is supported.\n");
313           return;
314      }
315
316      if (A_p->type & MATRIX_FORMAT_SYM) {
317        Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support symmetric storage scheme");
318        return;
319      }
320
321      /* open the file */
322      fileHandle_p = fopen(fileName_p, "w");
323      if (fileHandle_p==NULL) {
324        Paso_setError(IO_ERROR,"file could not be opened for writing");
325        return;
326      }
327
328      if (A_p->type & MATRIX_FORMAT_CSC) {
329        Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support CSC yet.");
330      } else {
331        mm_initialize_typecode(&matcode);
332        mm_set_matrix(&matcode);
333        mm_set_coordinate(&matcode);
334        mm_set_real(&matcode);
335
336        N= A_p->numRows;
337        M= A_p->numCols;
338        mm_write_banner(fileHandle_p, matcode);
339        mm_write_mtx_crd_size(fileHandle_p, N, M, A_p->pattern->ptr[N]);
340        if(A_p->row_block_size==1)
341          for (i=0; i<N; i++)
342            for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
343              fprintf(fileHandle_p, "%d %d %25.15e\n", i+1, A_p->pattern->index[iptr_ij]+1, A_p->val[iptr_ij]);
344            }
345
346        if(A_p->row_block_size==2)
347          for (i=0; i<N; i++) {
348            for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
349              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]);
350              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]);
351              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]);
352              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]);
353            }
354          }
355
356        if(A_p->row_block_size==3)
357          for (i=0; i<N; i++) {
358           for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
359              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]);
360              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]);
361              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]);
362              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]);
363              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]);
364              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]);
365              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]);
366              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]);
367              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]);
368            }
369          }
370
371      }
372
373      /* close the file */
374      fclose(fileHandle_p);
375
376      return;
377    }
378

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