# Contents of /trunk/paso/src/SparseMatrix.c

Revision 3312 - (show annotations)
Tue Oct 26 07:54:58 2010 UTC (8 years, 11 months ago) by gross
File MIME type: text/plain
File size: 12377 byte(s)
```last step for a clean up version of the AMG
```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2010 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 17 /* Paso: SparseMatrix */ 18 19 /**************************************************************/ 20 21 /* Author: Lutz Gross, l.gross@uq.edu.au */ 22 23 /**************************************************************/ 24 25 #include "Paso.h" 26 #include "SparseMatrix.h" 27 #include "MKL.h" 28 #include "Preconditioner.h" 29 #include "UMFPACK.h" 30 #include "TRILINOS.h" 31 #include "mmio.h" 32 33 /**************************************************************/ 34 static void swap( index_t*, index_t*, double*, int, int ); 35 static void q_sort( index_t*, index_t*, double*, int, int ); 36 /*static void print_entries( index_t*, index_t*, double* );*/ 37 38 static int M, N, nz; 39 40 41 /* debug: print the entries */ 42 /* 43 void print_entries( index_t *r, index_t *c, double *v ) 44 { 45 int i; 46 47 for( i=0; i begin ) 79 { 80 pivot = N * row[begin] + col[begin]; 81 l = begin + 1; 82 r = end; 83 84 while( l < r ) 85 { 86 lval = N * row[l] + col[l]; 87 if( lval < pivot ) 88 l++; 89 else 90 { 91 r--; 92 swap( row, col, val, l, r ); 93 } 94 } 95 l--; 96 swap( row, col, val, begin, l ); 97 q_sort( row, col, val, begin, l ); 98 q_sort( row, col, val, r, end ); 99 } 100 } 101 102 103 /* allocates a SparseMatrix of type type using the given matrix pattern 104 Values are initialized by zero. 105 if patternIsUnrolled and type & MATRIX_FORMAT_BLK1, it is assumed that the pattern is allready unrolled to match the requested block size 106 and offsets otherwise unrolling and offset adjustment will be performed. 107 */ 108 Paso_SparseMatrix* Paso_SparseMatrix_alloc(Paso_SparseMatrixType type,Paso_Pattern *pattern, int row_block_size, int col_block_size, const bool_t patternIsUnrolled) { 109 110 Paso_SparseMatrix*out=NULL; 111 Paso_SparseMatrixType pattern_format_out; 112 bool_t unroll=FALSE; 113 114 if (patternIsUnrolled) { 115 if (! XNOR(type & MATRIX_FORMAT_OFFSET1, pattern->type & MATRIX_FORMAT_OFFSET1) ) { 116 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_alloc: requested offset and pattern offset does not match."); 117 return NULL; 118 } 119 } 120 /* do we need to apply unrolling ? */ 121 unroll 122 /* we don't like non-square blocks */ 123 = (row_block_size!=col_block_size) 124 /* or any block size bigger than 3 */ 125 || (col_block_size>3) 126 /* or if lock size one requested and the block size is not 1 */ 127 || ((type & MATRIX_FORMAT_BLK1) && (col_block_size>1) ) 128 /* offsets don't match */ 129 || ( (type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & MATRIX_FORMAT_OFFSET1) ) ; 130 131 pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? MATRIX_FORMAT_OFFSET1: MATRIX_FORMAT_DEFAULT; 132 133 Esys_resetError(); 134 out=MEMALLOC(1,Paso_SparseMatrix); 135 if (! Esys_checkPtr(out)) { 136 out->pattern=NULL; 137 out->val=NULL; 138 out->reference_counter=1; 139 out->type=type; 140 out->solver_package=PASO_PASO; 141 out->solver_p=NULL; 142 143 /* ====== compressed sparse columns === */ 144 if (type & MATRIX_FORMAT_CSC) { 145 if (unroll) { 146 if (patternIsUnrolled) { 147 out->pattern=Paso_Pattern_getReference(pattern); 148 } else { 149 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size); 150 } 151 out->row_block_size=1; 152 out->col_block_size=1; 153 } else { 154 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1); 155 out->row_block_size=row_block_size; 156 out->col_block_size=col_block_size; 157 } 158 if (Esys_noError()) { 159 out->numRows = out->pattern->numInput; 160 out->numCols = out->pattern->numOutput; 161 } 162 163 } else { 164 /* ====== compressed sparse row === */ 165 if (unroll) { 166 if (patternIsUnrolled) { 167 out->pattern=Paso_Pattern_getReference(pattern); 168 } else { 169 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size); 170 } 171 out->row_block_size=1; 172 out->col_block_size=1; 173 } else { 174 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1); 175 out->row_block_size=row_block_size; 176 out->col_block_size=col_block_size; 177 } 178 if (Esys_noError()) { 179 out->numRows = out->pattern->numOutput; 180 out->numCols = out->pattern->numInput; 181 } 182 } 183 if (Esys_noError()) { 184 if (type & MATRIX_FORMAT_DIAGONAL_BLOCK) { 185 out->block_size=MIN(out->row_block_size,out->col_block_size); 186 } else { 187 out->block_size=out->row_block_size*out->col_block_size; 188 } 189 out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size); 190 191 out->val=MEMALLOC(out->len,double); 192 if (! Esys_checkPtr(out->val)) Paso_SparseMatrix_setValues(out,DBLE(0)); 193 } 194 } 195 /* all done: */ 196 if (Esys_noError()) { 197 return out; 198 } else { 199 Paso_SparseMatrix_free(out); 200 return NULL; 201 } 202 } 203 204 /* returns a reference to Paso_SparseMatrix in */ 205 206 Paso_SparseMatrix* Paso_SparseMatrix_getReference(Paso_SparseMatrix* in) { 207 if (in!=NULL) ++(in->reference_counter); 208 return in; 209 } 210 211 /* deallocates a SparseMatrix: */ 212 213 void Paso_SparseMatrix_free(Paso_SparseMatrix* in) { 214 if (in!=NULL) { 215 in->reference_counter--; 216 if (in->reference_counter<=0) { 217 switch(in->solver_package) { 218 219 case PASO_SMOOTHER: 220 Paso_Preconditioner_LocalSmoother_free((Paso_Preconditioner_LocalSmoother*) in->solver_p); 221 break; 222 223 case PASO_MKL: 224 Paso_MKL_free(in); /* releases solver_p */ 225 break; 226 227 case PASO_UMFPACK: 228 Paso_UMFPACK_free(in); /* releases solver_p */ 229 break; 230 } 231 MEMFREE(in->val); 232 Paso_Pattern_free(in->pattern); 233 MEMFREE(in); 234 } 235 } 236 } 237 238 Paso_SparseMatrix* Paso_SparseMatrix_loadMM_toCSR( char *fileName_p ) 239 { 240 index_t *col_ind = NULL; 241 index_t *row_ind = NULL; 242 index_t *row_ptr = NULL; 243 double *val = NULL; 244 FILE *fileHandle_p = NULL; 245 Paso_Pattern* mainPattern=NULL; 246 Paso_SparseMatrix *out = NULL; 247 int i, curr_row, scan_ret; 248 MM_typecode matrixCode; 249 Esys_resetError(); 250 251 /* open the file */ 252 fileHandle_p = fopen( fileName_p, "r" ); 253 if( fileHandle_p == NULL ) 254 { 255 Esys_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Cannot read file for reading."); 256 return NULL; 257 } 258 259 /* process banner */ 260 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 ) 261 { 262 Esys_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Error processing MM banner."); 263 fclose( fileHandle_p ); 264 return NULL; 265 } 266 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) ) 267 { 268 269 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_loadMM_toCSR: found Matrix Market type is not supported."); 270 fclose( fileHandle_p ); 271 return NULL; 272 } 273 274 /* get matrix size */ 275 if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 ) 276 { 277 Esys_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not parse matrix size"); 278 fclose( fileHandle_p ); 279 return NULL; 280 } 281 282 /* prepare storage */ 283 col_ind = MEMALLOC( nz, index_t ); 284 row_ind = MEMALLOC( nz, index_t ); 285 val = MEMALLOC( nz, double ); 286 287 row_ptr = MEMALLOC( (M+1), index_t ); 288 289 if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL ) 290 { 291 Esys_setError(MEMORY_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not allocate memory" ); 292 fclose( fileHandle_p ); 293 return NULL; 294 } 295 296 /* perform actual read of elements */ 297 for( i=0; ival[i] = val[i]; 331 332 Paso_Pattern_free(mainPattern); 333 MEMFREE( val ); 334 MEMFREE( row_ind ); 335 return out; 336 } 337 338 339 void Paso_SparseMatrix_saveMM(Paso_SparseMatrix * A_p, char * fileName_p) { 340 FILE * fileHandle_p = NULL; 341 dim_t N,M,i,j, irow, icol, ib, irb, icb, iptr; 342 MM_typecode matcode; 343 const dim_t col_block_size = A_p->col_block_size; 344 const dim_t row_block_size = A_p->row_block_size; 345 const dim_t block_size = A_p->block_size; 346 347 if (col_block_size !=row_block_size) { 348 Esys_setError(TYPE_ERROR, "Paso_SparseMatrix_saveMM: currently only square blocks are supported."); 349 return; 350 } 351 352 /* open the file */ 353 fileHandle_p = fopen(fileName_p, "w"); 354 if (fileHandle_p==NULL) { 355 Esys_setError(IO_ERROR,"file could not be opened for writing"); 356 return; 357 } 358 if (A_p->type & MATRIX_FORMAT_CSC) { 359 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support CSC yet."); 360 } else { 361 mm_initialize_typecode(&matcode); 362 mm_set_matrix(&matcode); 363 mm_set_coordinate(&matcode); 364 mm_set_real(&matcode); 365 366 N=Paso_SparseMatrix_getNumRows(A_p); 367 M=Paso_SparseMatrix_getNumCols(A_p); 368 mm_write_banner(fileHandle_p, matcode); 369 mm_write_mtx_crd_size(fileHandle_p, N*row_block_size, M*col_block_size, A_p->pattern->ptr[N]*block_size); 370 371 if (A_p->type & MATRIX_FORMAT_DIAGONAL_BLOCK) { 372 for (i=0; ipattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 374 j=A_p->pattern->index[iptr]; 375 for (ib=0;ibval[iptr*block_size+ib]); 379 } 380 } 381 } 382 } else { 383 for (i=0; ipattern->ptr[i];iptrpattern->ptr[i+1]; ++iptr) { 385 j=A_p->pattern->index[iptr]; 386 for (irb=0;irbval[iptr*block_size+irb+row_block_size*icb]); 391 } 392 } 393 } 394 } 395 } 396 } 397 /* close the file */ 398 fclose(fileHandle_p); 399 return; 400 } 401 402 index_t* Paso_SparseMatrix_borrowMainDiagonalPointer(Paso_SparseMatrix * A_p) 403 { 404 return Paso_Pattern_borrowMainDiagonalPointer(A_p->pattern); 405 } 406 407 dim_t Paso_SparseMatrix_getNumColors(Paso_SparseMatrix* A_p) 408 { 409 return Paso_Pattern_getNumColors(A_p->pattern); 410 } 411 index_t* Paso_SparseMatrix_borrowColoringPointer(Paso_SparseMatrix* A_p) 412 { 413 return Paso_Pattern_borrowColoringPointer(A_p->pattern); 414 } 415 dim_t Paso_SparseMatrix_maxDeg(Paso_SparseMatrix * A_p) 416 { 417 return Paso_Pattern_maxDeg(A_p->pattern); 418 } 419 dim_t Paso_SparseMatrix_getTotalNumRows(const Paso_SparseMatrix* A){ 420 return A->numRows * A->row_block_size; 421 } 422 423 dim_t Paso_SparseMatrix_getTotalNumCols(const Paso_SparseMatrix* A){ 424 return A->numCols * A->col_block_size; 425 } 426 dim_t Paso_SparseMatrix_getNumRows(Paso_SparseMatrix* A) { 427 return A->numRows; 428 } 429 dim_t Paso_SparseMatrix_getNumCols(Paso_SparseMatrix* A) { 430 return A->numCols; 431 }

 ViewVC Help Powered by ViewVC 1.1.26