1 #include "petsc.h" 2 #include "../src/mat/impls/aij/seq/aij.h" 3 #include "spbas.h" 4 5 6 7 8 /* 9 spbas_memory_requirement: 10 Calculate the number of bytes needed to store tha matrix 11 */ 12 #undef __FUNCT__ 13 #define __FUNCT__ "spbas_memory_requirement" 14 long int spbas_memory_requirement( spbas_matrix matrix) 15 { 16 17 18 // Requirement for 19 long int memreq = 6 * sizeof(PetscInt) + // nrows, ncols, nnz, n_alloc_icol, 20 // n_alloc_val, col_idx_type 21 sizeof(PetscTruth) + // block_data 22 sizeof(PetscScalar**) + // values 23 sizeof(PetscScalar*) + // alloc_val 24 2 * sizeof(int**) + // icols, icols0 25 2 * sizeof(PetscInt*) + // row_nnz, alloc_icol 26 matrix.nrows * sizeof(PetscInt) + // row_nnz[*] 27 matrix.nrows * sizeof(PetscInt*); // icols[*] 28 29 // icol0[*] 30 if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) 31 { memreq += matrix.nrows * sizeof(PetscInt); } 32 33 // icols[*][*] 34 if (matrix.block_data) 35 { memreq += matrix.n_alloc_icol * sizeof(PetscInt); } 36 else 37 { memreq += matrix.nnz * sizeof(PetscInt); } 38 39 if (matrix.values) 40 { 41 memreq += matrix.nrows * sizeof(PetscScalar*); // values[*] 42 // values[*][*] 43 if (matrix.block_data) 44 { memreq += matrix.n_alloc_val * sizeof(PetscScalar); } 45 else 46 { memreq += matrix.nnz * sizeof(PetscScalar); } 47 } 48 49 return memreq; 50 } 51 52 /* 53 spbas_allocate_pattern: 54 allocate the pattern arrays row_nnz, icols and optionally values 55 */ 56 #undef __FUNCT__ 57 #define __FUNCT__ "spbas_allocate_pattern" 58 PetscErrorCode spbas_allocate_pattern( spbas_matrix * result, PetscTruth do_values) 59 { 60 PetscErrorCode ierr; 61 PetscInt nrows = result->nrows; 62 PetscInt col_idx_type = result->col_idx_type; 63 64 PetscFunctionBegin; 65 // Allocate sparseness pattern 66 ierr = PetscMalloc(nrows*sizeof(PetscInt),&result->row_nnz); CHKERRQ(ierr); 67 ierr = PetscMalloc(nrows*sizeof(PetscInt*),&result->icols); CHKERRQ(ierr); 68 69 // If offsets are given wrt an array, create array 70 if (col_idx_type == SPBAS_OFFSET_ARRAY) 71 { 72 ierr = PetscMalloc(nrows*sizeof(PetscInt),&result->icol0); CHKERRQ(ierr); 73 } 74 else 75 { 76 result->icol0 = NULL; 77 } 78 79 // If values are given, allocate values array 80 if (do_values) 81 { 82 ierr = PetscMalloc(nrows*sizeof(PetscScalar*),&result->values); 83 CHKERRQ(ierr); 84 } 85 else 86 { 87 result->values = NULL; 88 } 89 90 PetscFunctionReturn(0); 91 } 92 93 /* 94 spbas_allocate_data: 95 in case of block_data: 96 Allocate the data arrays alloc_icol and optionally alloc_val, 97 set appropriate pointers from icols and values; 98 in case of !block_data: 99 Allocate the arrays icols[i] and optionally values[i] 100 */ 101 #undef __FUNCT__ 102 #define __FUNCT__ "spbas_allocate_data" 103 PetscErrorCode spbas_allocate_data( spbas_matrix * result) 104 { 105 PetscInt i; 106 PetscInt nnz = result->nnz; 107 PetscInt nrows = result->nrows; 108 PetscInt r_nnz; 109 PetscErrorCode ierr; 110 PetscTruth do_values = (result->values != NULL) ? PETSC_TRUE : PETSC_FALSE; 111 PetscTruth block_data = result->block_data; 112 113 PetscFunctionBegin; 114 if (block_data) 115 { 116 // Allocate the column number array and point to it 117 result->n_alloc_icol = nnz; 118 ierr = PetscMalloc(nnz*sizeof(PetscInt), &result->alloc_icol); 119 CHKERRQ(ierr); 120 result->icols[0]= result->alloc_icol; 121 for (i=1; i<nrows; i++) 122 { 123 result->icols[i] = result->icols[i-1] + result->row_nnz[i-1]; 124 } 125 126 // Allocate the value array and point to it 127 if (do_values) 128 { 129 result->n_alloc_val= nnz; 130 ierr = PetscMalloc(nnz*sizeof(PetscScalar), &result->alloc_val); 131 CHKERRQ(ierr); 132 result->values[0]= result->alloc_val; 133 for (i=1; i<nrows; i++) 134 { 135 result->values[i] = result->values[i-1] + result->row_nnz[i-1]; 136 } 137 } 138 } 139 else 140 { 141 for (i=0; i<nrows; i++) 142 { 143 r_nnz = result->row_nnz[i]; 144 ierr = PetscMalloc(r_nnz*sizeof(PetscInt), &result->icols[i]); 145 CHKERRQ(ierr); 146 } 147 if (do_values) 148 { 149 for (i=0; i<nrows; i++) 150 { 151 r_nnz = result->row_nnz[i]; 152 ierr = PetscMalloc(r_nnz*sizeof(PetscScalar), 153 &result->values[i]); CHKERRQ(ierr); 154 } 155 } 156 } 157 158 PetscFunctionReturn(0); 159 } 160 161 /* 162 spbas_row_order_icol 163 determine if row i1 should come 164 + before row i2 in the sorted rows (return -1), 165 + after (return 1) 166 + is identical (return 0). 167 */ 168 #undef __FUNCT__ 169 #define __FUNCT__ "spbas_row_order_icol" 170 int spbas_row_order_icol(PetscInt i1, PetscInt i2, PetscInt *irow_in, PetscInt *icol_in,PetscInt col_idx_type) 171 { 172 PetscInt j; 173 PetscInt nnz1 = irow_in[i1+1] - irow_in[i1]; 174 PetscInt nnz2 = irow_in[i2+1] - irow_in[i2]; 175 176 PetscInt * icol1 = &icol_in[irow_in[i1]]; 177 PetscInt * icol2 = &icol_in[irow_in[i2]]; 178 179 if (nnz1<nnz2) {return -1;} 180 if (nnz1>nnz2) {return 1;} 181 182 if (col_idx_type == SPBAS_COLUMN_NUMBERS) 183 { 184 for (j=0; j<nnz1; j++) 185 { 186 if (icol1[j]< icol2[j]) {return -1;} 187 if (icol1[j]> icol2[j]) {return 1;} 188 } 189 } 190 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 191 { 192 for (j=0; j<nnz1; j++) 193 { 194 if (icol1[j]-i1< icol2[j]-i2) {return -1;} 195 if (icol1[j]-i1> icol2[j]-i2) {return 1;} 196 } 197 } 198 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 199 { 200 for (j=1; j<nnz1; j++) 201 { 202 if (icol1[j]-icol1[0] < icol2[j]-icol2[0]) {return -1;} 203 if (icol1[j]-icol1[0] > icol2[j]-icol2[0]) {return 1;} 204 } 205 } 206 return 0; 207 } 208 209 /* 210 spbas_mergesort_icols: 211 return a sorting of the rows in which identical sparseness patterns are 212 next to each other 213 */ 214 #undef __FUNCT__ 215 #define __FUNCT__ "spbas_mergesort_icols" 216 PetscErrorCode spbas_mergesort_icols( PetscInt nrows, PetscInt * irow_in, PetscInt * icol_in,PetscInt col_idx_type, PetscInt *isort) 217 { 218 PetscErrorCode ierr; 219 PetscInt istep; // Chunk-sizes of already sorted parts of arrays 220 221 PetscInt i, i1, i2; // Loop counters for (partly) sorted arrays 222 223 PetscInt istart, i1end, i2end; // start of newly sorted array part, end of both 224 // parts 225 226 PetscInt *ialloc; // Allocated arrays 227 PetscInt *iswap; // auxiliary pointers for swapping 228 PetscInt *ihlp1; // Pointers to new version of arrays, 229 PetscInt *ihlp2; // Pointers to previous version of arrays, 230 231 PetscFunctionBegin; 232 233 // Create arrays 234 ierr = PetscMalloc(nrows*sizeof(PetscInt),&ialloc); CHKERRQ(ierr); 235 236 ihlp1 = ialloc; 237 ihlp2 = isort; 238 239 // Sorted array chunks are first 1 long, and increase until they are the 240 // complete array 241 for (istep=1; istep<nrows; istep*=2) 242 { 243 // 244 // Combine sorted parts 245 // istart:istart+istep-1 and istart+istep-1:istart+2*istep-1 246 // of ihlp2 and vhlp2 247 // 248 // into one sorted part 249 // istart:istart+2*istep-1 250 // of ihlp1 and vhlp1 251 // 252 for (istart=0; istart<nrows; istart+=2*istep) 253 { 254 255 // Set counters and bound array part endings 256 i1=istart; i1end = i1+istep; if (i1end>nrows) {i1end=nrows;} 257 i2=istart+istep; i2end = i2+istep; if (i2end>nrows) {i2end=nrows;} 258 259 // Merge the two array parts 260 for (i=istart; i<i2end; i++) 261 { 262 if (i1<i1end && i2<i2end && 263 spbas_row_order_icol( ihlp2[i1], ihlp2[i2], 264 irow_in, icol_in, col_idx_type) < 0) 265 { 266 ihlp1[i] = ihlp2[i1]; 267 i1++; 268 } 269 else if (i2<i2end ) 270 { 271 ihlp1[i] = ihlp2[i2]; 272 i2++; 273 } 274 else 275 { 276 ihlp1[i] = ihlp2[i1]; 277 i1++; 278 } 279 } 280 } 281 282 // Swap the two array sets 283 iswap = ihlp2; ihlp2 = ihlp1; ihlp1 = iswap; 284 } 285 286 // Copy one more time in case the sorted arrays are the temporary ones 287 if (ihlp2 != isort) 288 { 289 for (i=0; i<nrows; i++) { isort[i] = ihlp2[i]; } 290 } 291 292 // Remove help arrays 293 ierr = PetscFree(ialloc); CHKERRQ(ierr); 294 PetscFunctionReturn(0); 295 } 296 297 298 299 /* 300 spbas_compress_pattern: 301 calculate a compressed sparseness pattern for a sparseness pattern 302 given in compressed row storage. The compressed sparseness pattern may 303 require (much) less memory. 304 */ 305 #undef __FUNCT__ 306 #define __FUNCT__ "spbas_compress_pattern" 307 PetscErrorCode spbas_compress_pattern(PetscInt *irow_in, PetscInt *icol_in, PetscInt nrows, PetscInt ncols, 308 PetscInt col_idx_type, spbas_matrix *B,PetscScalar *mem_reduction) 309 { 310 PetscInt nnz = irow_in[nrows]; 311 long int mem_orig = (nrows + nnz) * sizeof(PetscInt); 312 long int mem_compressed; 313 PetscErrorCode ierr; 314 PetscInt *isort; 315 PetscInt *icols; 316 PetscInt row_nnz; 317 PetscInt *ipoint; 318 PetscTruth *used; 319 PetscInt ptr; 320 PetscInt i,j; 321 322 const PetscTruth no_values = PETSC_FALSE; 323 324 PetscFunctionBegin; 325 326 // Allocate the structure of the new matrix 327 B->nrows = nrows; 328 B->ncols = ncols; 329 B->nnz = nnz; 330 B->col_idx_type= col_idx_type; 331 B->block_data = PETSC_TRUE; 332 ierr = spbas_allocate_pattern( B, no_values); CHKERRQ(ierr); 333 334 // When using an offset array, set it 335 if (col_idx_type==SPBAS_OFFSET_ARRAY) 336 { 337 for (i=0; i<nrows; i++) {B->icol0[i] = icol_in[irow_in[i]];} 338 } 339 340 // Allocate the ordering for the rows 341 ierr = PetscMalloc(nrows*sizeof(PetscInt),&isort); CHKERRQ(ierr); 342 ierr = PetscMalloc(nrows*sizeof(PetscInt),&ipoint); CHKERRQ(ierr); 343 ierr = PetscMalloc(nrows*sizeof(PetscTruth),&used); CHKERRQ(ierr); 344 345 // Initialize the sorting 346 memset((void*) used, 0, nrows*sizeof(PetscTruth)); 347 for (i = 0; i<nrows; i++) 348 { 349 B->row_nnz[i] = irow_in[i+1]-irow_in[i]; 350 isort[i] = i; 351 ipoint[i]= i; 352 } 353 354 // Sort the rows so that identical columns will be next to each other 355 ierr = spbas_mergesort_icols( nrows, irow_in, icol_in, col_idx_type, isort); CHKERRQ(ierr); 356 printf("Rows have been sorted for patterns\n"); 357 358 // Replace identical rows with the first one in the list 359 for (i=1; i<nrows; i++) 360 { 361 if (spbas_row_order_icol(isort[i-1], isort[i], 362 irow_in, icol_in, col_idx_type) == 0) 363 { 364 ipoint[isort[i]] = ipoint[isort[i-1]]; 365 } 366 } 367 368 // Collect the rows which are used 369 for (i=0; i<nrows; i++) {used[ipoint[i]] = PETSC_TRUE;} 370 371 // Calculate needed memory 372 B->n_alloc_icol = 0; 373 for (i=0; i<nrows; i++) 374 { 375 if (used[i]) {B->n_alloc_icol += B->row_nnz[i];} 376 } 377 378 ierr = PetscMalloc(B->n_alloc_icol*sizeof(PetscInt),&B->alloc_icol); 379 CHKERRQ(ierr); 380 381 // Fill in the diagonal offsets for the rows which store their own data 382 ptr = 0; 383 for (i=0; i<B->nrows; i++) 384 { 385 if (used[i]) 386 { 387 B->icols[i] = &B->alloc_icol[ptr]; 388 icols = &icol_in[irow_in[i]]; 389 row_nnz = B->row_nnz[i]; 390 if (col_idx_type == SPBAS_COLUMN_NUMBERS) 391 { 392 for (j=0; j<row_nnz; j++) 393 { 394 B->icols[i][j] = icols[j]; 395 } 396 } 397 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 398 { 399 for (j=0; j<row_nnz; j++) 400 { 401 B->icols[i][j] = icols[j]-i; 402 } 403 } 404 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 405 { 406 for (j=0; j<row_nnz; j++) 407 { 408 B->icols[i][j] = icols[j]-icols[0]; 409 } 410 } 411 ptr += B->row_nnz[i]; 412 } 413 } 414 415 // Point to the right places for all data 416 for (i=0; i<nrows; i++) 417 { 418 B->icols[i] = B->icols[ipoint[i]]; 419 } 420 printf("Row patterns have been compressed\n"); 421 printf(" (%6.2f nonzeros per row)\n", 422 (PetscScalar) nnz / (PetscScalar) nrows); 423 424 425 // Clean up 426 ierr=PetscFree(isort); CHKERRQ(ierr); 427 ierr=PetscFree(used); CHKERRQ(ierr); 428 ierr=PetscFree(ipoint); CHKERRQ(ierr); 429 430 mem_compressed = spbas_memory_requirement( *B ); 431 *mem_reduction = 100.0 * (PetscScalar)(mem_orig-mem_compressed)/ 432 (PetscScalar) mem_orig; 433 PetscFunctionReturn(0); 434 } 435 436 /* 437 spbas_incomplete_cholesky 438 Incomplete Cholesky decomposition 439 */ 440 #include "spbas_cholesky.h" 441 442 /* 443 spbas_delete : de-allocate the arrays owned by this matrix 444 */ 445 #undef __FUNCT__ 446 #define __FUNCT__ "spbas_delete" 447 PetscErrorCode spbas_delete(spbas_matrix matrix) 448 { 449 PetscInt i; 450 PetscErrorCode ierr; 451 PetscFunctionBegin; 452 if (matrix.block_data) 453 { 454 ierr=PetscFree(matrix.alloc_icol); CHKERRQ(ierr) 455 if (matrix.values){ierr=PetscFree(matrix.alloc_val);CHKERRQ(ierr);} 456 } 457 else 458 { 459 for (i=0; i<matrix.nrows; i++) 460 { ierr=PetscFree(matrix.icols[i]); CHKERRQ(ierr);} 461 ierr = PetscFree(matrix.icols);CHKERRQ(ierr); 462 if (matrix.values) 463 { 464 for (i=0; i<matrix.nrows; i++) 465 { ierr=PetscFree(matrix.values[i]); CHKERRQ(ierr);} 466 } 467 } 468 469 ierr=PetscFree(matrix.row_nnz); CHKERRQ(ierr); 470 ierr=PetscFree(matrix.icols); CHKERRQ(ierr); 471 if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) 472 {ierr=PetscFree(matrix.icol0);CHKERRQ(ierr);} 473 if (matrix.values) 474 {ierr=PetscFree(matrix.values);CHKERRQ(ierr);} 475 476 PetscFunctionReturn(0); 477 } 478 479 /* 480 spbas_matrix_to_crs: 481 Convert an spbas_matrix to compessed row storage 482 */ 483 #undef __FUNCT__ 484 #define __FUNCT__ "spbas_matrix_to_crs" 485 PetscErrorCode spbas_matrix_to_crs(spbas_matrix matrix_A,MatScalar **val_out, PetscInt **irow_out, PetscInt **icol_out) 486 { 487 PetscInt nrows = matrix_A.nrows; 488 PetscInt nnz = matrix_A.nnz; 489 PetscInt i,j,r_nnz,i0; 490 PetscInt *irow; 491 PetscInt *icol; 492 PetscInt *icol_A; 493 MatScalar *val; 494 PetscScalar *val_A; 495 PetscInt col_idx_type = matrix_A.col_idx_type; 496 PetscTruth do_values = matrix_A.values != NULL; 497 PetscErrorCode ierr; 498 499 PetscFunctionBegin; 500 ierr = PetscMalloc( sizeof(PetscInt) * (nrows+1), &irow); CHKERRQ(ierr); 501 ierr = PetscMalloc( sizeof(PetscInt) * nnz, &icol); CHKERRQ(ierr); 502 *icol_out = icol; *irow_out=irow; 503 if (do_values) 504 { 505 ierr = PetscMalloc( sizeof(MatScalar) * nnz, &val); CHKERRQ(ierr); 506 *val_out = val; *icol_out = icol; *irow_out=irow; 507 } 508 509 irow[0]=0; 510 for (i=0; i<nrows; i++) 511 { 512 r_nnz = matrix_A.row_nnz[i]; 513 i0 = irow[i]; 514 irow[i+1] = i0 + r_nnz; 515 icol_A = matrix_A.icols[i]; 516 517 if (do_values) 518 { 519 val_A = matrix_A.values[i]; 520 for (j=0; j<r_nnz; j++) 521 { 522 icol[i0+j] = icol_A[j]; 523 val[i0+j] = val_A[j]; 524 } 525 } 526 else 527 { 528 for (j=0; j<r_nnz; j++) { icol[i0+j] = icol_A[j]; } 529 } 530 531 if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 532 { 533 for (j=0; j<r_nnz; j++) { icol[i0+j] += i; } 534 } 535 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 536 { 537 i0 = matrix_A.icol0[i]; 538 for (j=0; j<r_nnz; j++) { icol[i0+j] += i0; } 539 } 540 541 } 542 543 PetscFunctionReturn(0); 544 } 545 546 /* 547 spbas_dump 548 Print the matrix in i,j,val-format 549 */ 550 #undef __FUNCT__ 551 #define __FUNCT__ "spbas_dump" 552 PetscErrorCode spbas_dump( const char *fname, spbas_matrix in_matrix) 553 { 554 FILE *outfile = fopen(fname, "w"); 555 PetscInt col_idx_type = in_matrix.col_idx_type; 556 PetscInt nrows=in_matrix.nrows; 557 PetscInt *icol; 558 PetscScalar *val; 559 PetscInt r_nnz, icol0; 560 PetscInt i,j; 561 PetscInt nwrite=1; 562 563 PetscFunctionBegin; 564 if (!outfile) 565 { 566 SETERRQ1( PETSC_ERR_FILE_OPEN, 567 "Error in spbas_dump: cannot open file '%s'\n",fname); 568 } 569 if (in_matrix.values) 570 { 571 for (i=0; i<nrows; i++) 572 { 573 icol = in_matrix.icols[i]; 574 val = in_matrix.values[i]; 575 r_nnz= in_matrix.row_nnz[i]; 576 if (col_idx_type == SPBAS_COLUMN_NUMBERS) 577 { 578 for (j=0; nwrite>0 && j<r_nnz; j++) 579 { nwrite = fprintf(outfile,"%d %d %e\n",i,icol[j],val[j]); } 580 } 581 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 582 { 583 for (j=0; nwrite>0 && j<r_nnz; j++) 584 { nwrite = fprintf(outfile,"%d %d %e\n",i,i+icol[j],val[j]);} 585 } 586 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 587 { 588 icol0 = in_matrix.icol0[i]; 589 for (j=0; nwrite>0 && j<r_nnz; j++) 590 { nwrite = fprintf(outfile,"%d %d %e\n",i,icol0+icol[j],val[j]); } 591 } 592 593 } 594 } 595 else 596 { 597 for (i=0; i<nrows; i++) 598 { 599 icol = in_matrix.icols[i]; 600 r_nnz= in_matrix.row_nnz[i]; 601 nwrite = 2; 602 if (col_idx_type == SPBAS_COLUMN_NUMBERS) 603 { 604 for (j=0; nwrite>0 && j<r_nnz; j++) 605 { nwrite = fprintf(outfile,"%d %d\n",i,icol[j]); } 606 } 607 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 608 { 609 for (j=0; nwrite>0 && j<r_nnz; j++) 610 { nwrite = fprintf(outfile,"%d %d\n",i,i+icol[j]); } 611 } 612 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 613 { 614 icol0 = in_matrix.icol0[i]; 615 for (j=0; nwrite>0 && j<r_nnz; j++) 616 { nwrite = fprintf(outfile,"%d %d\n",i,icol0+icol[j]); } 617 } 618 } 619 } 620 if (nwrite<0) 621 { 622 SETERRQ1(PETSC_ERR_FILE_WRITE, 623 "Error in spbas_dump: cannot write to file '%s'\n",fname); 624 } 625 fclose(outfile); 626 PetscFunctionReturn(0); 627 } 628 629 /* 630 spbas_transpose 631 return the transpose of a matrix 632 */ 633 #undef __FUNCT__ 634 #define __FUNCT__ "spbas_transpose" 635 PetscErrorCode spbas_transpose( spbas_matrix in_matrix, spbas_matrix * result) 636 { 637 PetscInt col_idx_type = in_matrix.col_idx_type; 638 PetscInt nnz = in_matrix.nnz; 639 PetscInt ncols = in_matrix.nrows; 640 PetscInt nrows = in_matrix.ncols; 641 PetscInt i,j,k; 642 PetscInt r_nnz; 643 PetscInt *irow; 644 PetscInt icol0; 645 PetscScalar * val; 646 PetscErrorCode ierr; 647 648 PetscFunctionBegin; 649 650 // Copy input values 651 result->nrows = nrows; 652 result->ncols = ncols; 653 result->nnz = nnz; 654 result->col_idx_type = SPBAS_COLUMN_NUMBERS; 655 result->block_data = PETSC_TRUE; 656 657 // Allocate sparseness pattern 658 ierr = spbas_allocate_pattern(result, in_matrix.values != NULL); 659 CHKERRQ(ierr); 660 661 // Count the number of nonzeros in each row 662 for (i = 0; i<nrows; i++) { result->row_nnz[i] = 0; } 663 664 for (i=0; i<ncols; i++) 665 { 666 r_nnz = in_matrix.row_nnz[i]; 667 irow = in_matrix.icols[i]; 668 if (col_idx_type == SPBAS_COLUMN_NUMBERS) 669 { 670 for (j=0; j<r_nnz; j++) { result->row_nnz[irow[j]]++; } 671 } 672 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 673 { 674 for (j=0; j<r_nnz; j++) { result->row_nnz[i+irow[j]]++; } 675 } 676 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 677 { 678 icol0=in_matrix.icol0[i]; 679 for (j=0; j<r_nnz; j++) { result->row_nnz[icol0+irow[j]]++; } 680 } 681 } 682 683 // Set the pointers to the data 684 ierr = spbas_allocate_data(result); CHKERRQ(ierr); 685 686 // Reset the number of nonzeros in each row 687 for (i = 0; i<nrows; i++) { result->row_nnz[i] = 0; } 688 689 // Fill the data arrays 690 if (in_matrix.values) 691 { 692 for (i=0; i<ncols; i++) 693 { 694 r_nnz = in_matrix.row_nnz[i]; 695 irow = in_matrix.icols[i]; 696 val = in_matrix.values[i]; 697 698 if (col_idx_type == SPBAS_COLUMN_NUMBERS) {icol0=0;} 699 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {icol0=i;} 700 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 701 {icol0=in_matrix.icol0[i];} 702 for (j=0; j<r_nnz; j++) 703 { 704 k = icol0 + irow[j]; 705 result->icols[k][result->row_nnz[k]] = i; 706 result->values[k][result->row_nnz[k]] = val[j]; 707 result->row_nnz[k]++; 708 } 709 } 710 } 711 else 712 { 713 for (i=0; i<ncols; i++) 714 { 715 r_nnz = in_matrix.row_nnz[i]; 716 irow = in_matrix.icols[i]; 717 718 if (col_idx_type == SPBAS_COLUMN_NUMBERS) {icol0=0;} 719 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {icol0=i;} 720 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 721 {icol0=in_matrix.icol0[i];} 722 723 for (j=0; j<r_nnz; j++) 724 { 725 k = icol0 + irow[j]; 726 result->icols[k][result->row_nnz[k]] = i; 727 result->row_nnz[k]++; 728 } 729 } 730 } 731 732 // Set return value 733 PetscFunctionReturn(0); 734 } 735 736 /* 737 spbas_mergesort 738 739 mergesort for an array of intergers and an array of associated 740 reals 741 742 on output, icol[0..nnz-1] is increasing; 743 val[0..nnz-1] has undergone the same permutation as icol 744 745 NB: val may be NULL: in that case, only the integers are sorted 746 747 */ 748 #undef __FUNCT__ 749 #define __FUNCT__ "spbas_mergesort" 750 PetscErrorCode spbas_mergesort(PetscInt nnz, PetscInt *icol, PetscScalar *val) 751 { 752 PetscInt istep; // Chunk-sizes of already sorted parts of arrays 753 PetscInt i, i1, i2; // Loop counters for (partly) sorted arrays 754 PetscInt istart, i1end, i2end; // start of newly sorted array part, end of both parts 755 PetscInt *ialloc; // Allocated arrays 756 PetscScalar *valloc=NULL; 757 PetscInt *iswap; // auxiliary pointers for swapping 758 PetscScalar *vswap; 759 PetscInt *ihlp1; // Pointers to new version of arrays, 760 PetscScalar *vhlp1=NULL; // (arrays under construction) 761 PetscInt *ihlp2; // Pointers to previous version of arrays, 762 PetscScalar *vhlp2=NULL; 763 PetscErrorCode ierr; 764 765 // Create arrays 766 ierr = PetscMalloc(nnz*sizeof(PetscInt),&ialloc); CHKERRQ(ierr); 767 ihlp1 = ialloc; 768 ihlp2 = icol; 769 770 if (val) 771 { 772 ierr = PetscMalloc(nnz*sizeof(PetscScalar),&valloc); CHKERRQ(ierr); 773 vhlp1 = valloc; 774 vhlp2 = val; 775 } 776 777 778 // Sorted array chunks are first 1 long, and increase until they are the 779 // complete array 780 for (istep=1; istep<nnz; istep*=2) 781 { 782 // 783 // Combine sorted parts 784 // istart:istart+istep-1 and istart+istep-1:istart+2*istep-1 785 // of ihlp2 and vhlp2 786 // 787 // into one sorted part 788 // istart:istart+2*istep-1 789 // of ihlp1 and vhlp1 790 // 791 for (istart=0; istart<nnz; istart+=2*istep) 792 { 793 794 // Set counters and bound array part endings 795 i1=istart; i1end = i1+istep; if (i1end>nnz) {i1end=nnz;} 796 i2=istart+istep; i2end = i2+istep; if (i2end>nnz) {i2end=nnz;} 797 798 // Merge the two array parts 799 if (val) 800 { 801 for (i=istart; i<i2end; i++) 802 { 803 if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) 804 { 805 ihlp1[i] = ihlp2[i1]; 806 vhlp1[i] = vhlp2[i1]; 807 i1++; 808 } 809 else if (i2<i2end ) 810 { 811 ihlp1[i] = ihlp2[i2]; 812 vhlp1[i] = vhlp2[i2]; 813 i2++; 814 } 815 else 816 { 817 ihlp1[i] = ihlp2[i1]; 818 vhlp1[i] = vhlp2[i1]; 819 i1++; 820 } 821 } 822 } 823 else 824 { 825 for (i=istart; i<i2end; i++) 826 { 827 if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) 828 { 829 ihlp1[i] = ihlp2[i1]; 830 i1++; 831 } 832 else if (i2<i2end ) 833 { 834 ihlp1[i] = ihlp2[i2]; 835 i2++; 836 } 837 else 838 { 839 ihlp1[i] = ihlp2[i1]; 840 i1++; 841 } 842 } 843 } 844 } 845 846 // Swap the two array sets 847 iswap = ihlp2; ihlp2 = ihlp1; ihlp1 = iswap; 848 vswap = vhlp2; vhlp2 = vhlp1; vhlp1 = vswap; 849 } 850 851 // Copy one more time in case the sorted arrays are the temporary ones 852 if (ihlp2 != icol) 853 { 854 for (i=0; i<nnz; i++) { icol[i] = ihlp2[i]; } 855 if (val) 856 { 857 for (i=0; i<nnz; i++) { val[i] = vhlp2[i]; } 858 } 859 } 860 861 // Remove help arrays 862 ierr = PetscFree(ialloc); CHKERRQ(ierr); 863 if(val){ierr = PetscFree(valloc); CHKERRQ(ierr);} 864 PetscFunctionReturn(0); 865 } 866 867 /* 868 spbas_apply_reordering_rows: 869 apply the given reordering to the rows: matrix_A = matrix_A(perm,:); 870 */ 871 #undef __FUNCT__ 872 #define __FUNCT__ "spbas_apply_reordering_rows" 873 PetscErrorCode spbas_apply_reordering_rows(spbas_matrix *matrix_A, const PetscInt *permutation) 874 { 875 PetscInt i,j,ip; 876 PetscInt nrows=matrix_A->nrows; 877 PetscInt * row_nnz; 878 PetscInt **icols; 879 PetscTruth do_values = matrix_A->values != NULL; 880 PetscScalar **vals=NULL; 881 PetscErrorCode ierr; 882 883 PetscFunctionBegin; 884 if (matrix_A->col_idx_type != SPBAS_DIAGONAL_OFFSETS) 885 { 886 SETERRQ( PETSC_ERR_SUP_SYS, 887 "must have diagonal offsets in pattern\n"); 888 } 889 890 if (do_values) 891 { 892 ierr = PetscMalloc( sizeof(PetscScalar*)*nrows, &vals); CHKERRQ(ierr); 893 } 894 ierr = PetscMalloc( sizeof(PetscInt)*nrows, &row_nnz); CHKERRQ(ierr); 895 ierr = PetscMalloc( sizeof(PetscInt*)*nrows, &icols); CHKERRQ(ierr); 896 897 for (i=0; i<nrows;i++) 898 { 899 ip = permutation[i]; 900 if (do_values) {vals[i] = matrix_A->values[ip];} 901 icols[i] = matrix_A->icols[ip]; 902 row_nnz[i] = matrix_A->row_nnz[ip]; 903 for (j=0; j<row_nnz[i]; j++) { icols[i][j] += ip-i; } 904 } 905 906 if (do_values){ ierr = PetscFree(matrix_A->values); CHKERRQ(ierr);} 907 ierr = PetscFree(matrix_A->icols); CHKERRQ(ierr); 908 ierr = PetscFree(matrix_A->row_nnz); CHKERRQ(ierr); 909 910 if (do_values) { matrix_A->values = vals; } 911 matrix_A->icols = icols; 912 matrix_A->row_nnz = row_nnz; 913 914 PetscFunctionReturn(0); 915 } 916 917 918 /* 919 spbas_apply_reordering_cols: 920 apply the given reordering to the columns: matrix_A(:,perm) = matrix_A; 921 */ 922 #undef __FUNCT__ 923 #define __FUNCT__ "spbas_apply_reordering_cols" 924 PetscErrorCode spbas_apply_reordering_cols( spbas_matrix *matrix_A,const PetscInt *permutation) 925 { 926 PetscInt i,j; 927 PetscInt nrows=matrix_A->nrows; 928 PetscInt row_nnz; 929 PetscInt *icols; 930 PetscTruth do_values = matrix_A->values != NULL; 931 PetscScalar *vals=NULL; 932 933 PetscErrorCode ierr; 934 935 PetscFunctionBegin; 936 937 if (matrix_A->col_idx_type != SPBAS_DIAGONAL_OFFSETS) 938 { 939 SETERRQ( PETSC_ERR_SUP_SYS, 940 "must have diagonal offsets in pattern\n"); 941 } 942 943 for (i=0; i<nrows;i++) 944 { 945 icols = matrix_A->icols[i]; 946 row_nnz = matrix_A->row_nnz[i]; 947 if (do_values) { vals = matrix_A->values[i]; } 948 949 for (j=0; j<row_nnz; j++) 950 { 951 icols[j] = permutation[i+icols[j]]-i; 952 } 953 ierr = spbas_mergesort(row_nnz, icols, vals); CHKERRQ(ierr); 954 } 955 956 PetscFunctionReturn(0); 957 } 958 959 /* 960 spbas_apply_reordering: 961 apply the given reordering: matrix_A(perm,perm) = matrix_A; 962 */ 963 #undef __FUNCT__ 964 #define __FUNCT__ "spbas_apply_reordering" 965 PetscErrorCode spbas_apply_reordering( spbas_matrix *matrix_A, const PetscInt *permutation, const PetscInt * inv_perm) 966 { 967 PetscErrorCode ierr; 968 PetscFunctionBegin; 969 ierr = spbas_apply_reordering_rows( matrix_A, inv_perm); CHKERRQ(ierr); 970 ierr = spbas_apply_reordering_cols( matrix_A, permutation); CHKERRQ(ierr); 971 PetscFunctionReturn(0); 972 } 973 974 #undef __FUNCT__ 975 #define __FUNCT__ "spbas_pattern_only" 976 PetscErrorCode spbas_pattern_only(PetscInt nrows, PetscInt ncols, PetscInt *ai, PetscInt *aj, spbas_matrix * result) 977 { 978 spbas_matrix retval; 979 PetscInt i, j, i0, r_nnz; 980 PetscErrorCode ierr; 981 982 PetscFunctionBegin; 983 984 // Copy input values 985 retval.nrows = nrows; 986 retval.ncols = ncols; 987 retval.nnz = ai[nrows]; 988 989 retval.block_data = PETSC_TRUE; 990 retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS; 991 992 // Allocate output matrix 993 ierr = spbas_allocate_pattern(&retval, PETSC_FALSE); CHKERRQ(ierr); 994 for (i=0; i<nrows; i++) {retval.row_nnz[i] = ai[i+1]-ai[i];} 995 ierr = spbas_allocate_data(&retval); CHKERRQ(ierr); 996 // Copy the structure 997 for (i = 0; i<retval.nrows; i++) 998 { 999 i0 = ai[i]; 1000 r_nnz = ai[i+1]-i0; 1001 1002 for (j=0; j<r_nnz; j++) 1003 { 1004 retval.icols[i][j] = aj[i0+j]-i; 1005 } 1006 } 1007 1008 // Set return value 1009 *result = retval; 1010 PetscFunctionReturn(0); 1011 } 1012 1013 1014 /* 1015 spbas_mark_row_power: 1016 Mark the columns in row 'row' which are nonzero in 1017 matrix^2log(marker). 1018 */ 1019 #undef __FUNCT__ 1020 #define __FUNCT__ "spbas_mark_row_power" 1021 PetscErrorCode spbas_mark_row_power( 1022 PetscInt *iwork, // marker-vector 1023 PetscInt row, // row for which the columns are marked 1024 spbas_matrix * in_matrix, // matrix for which the power is being 1025 // calculated 1026 PetscInt marker, // marker-value: 2^power 1027 PetscInt minmrk, // lower bound for marked points 1028 PetscInt maxmrk) // upper bound for marked points 1029 { 1030 PetscErrorCode ierr; 1031 PetscInt i,j, nnz; 1032 1033 PetscFunctionBegin; 1034 nnz = in_matrix->row_nnz[row]; 1035 1036 // For higher powers, call this function recursively 1037 if (marker>1) 1038 { 1039 for (i=0; i<nnz;i++) 1040 { 1041 j = row + in_matrix->icols[row][i]; 1042 if (minmrk<=j && j<maxmrk && iwork[j] < marker ) 1043 { 1044 ierr = spbas_mark_row_power( iwork, row + in_matrix->icols[row][i], 1045 in_matrix, marker/2,minmrk,maxmrk); 1046 CHKERRQ(ierr); 1047 iwork[j] |= marker; 1048 } 1049 } 1050 } 1051 else 1052 { 1053 // Mark the columns reached 1054 for (i=0; i<nnz;i++) 1055 { 1056 j = row + in_matrix->icols[row][i]; 1057 if (minmrk<=j && j<maxmrk) { iwork[j] |= 1; } 1058 } 1059 } 1060 1061 PetscFunctionReturn(0); 1062 } 1063 1064 1065 /* 1066 spbas_power 1067 Calculate sparseness patterns for incomplete Cholesky decompositions 1068 of a given order: (almost) all nonzeros of the matrix^(order+1) which 1069 are inside the band width are found and stored in the output sparseness 1070 pattern. 1071 */ 1072 #undef __FUNCT__ 1073 #define __FUNCT__ "spbas_power" 1074 PetscErrorCode spbas_power (spbas_matrix in_matrix,PetscInt power, spbas_matrix * result) 1075 { 1076 spbas_matrix retval; 1077 PetscInt nrows = in_matrix.nrows; 1078 PetscInt ncols = in_matrix.ncols; 1079 PetscInt i, j, kend; 1080 PetscInt nnz, inz; 1081 PetscInt *iwork; 1082 PetscInt marker; 1083 PetscInt maxmrk=0; 1084 PetscErrorCode ierr; 1085 1086 PetscFunctionBegin; 1087 if (in_matrix.col_idx_type != SPBAS_DIAGONAL_OFFSETS) 1088 { 1089 SETERRQ( PETSC_ERR_SUP_SYS, 1090 "must have diagonal offsets in pattern\n"); 1091 } 1092 1093 if (ncols != nrows) 1094 { 1095 SETERRQ( PETSC_ERR_ARG_INCOMP, 1096 "Dimension error\n"); 1097 } 1098 1099 if (in_matrix.values) 1100 { 1101 SETERRQ( PETSC_ERR_ARG_INCOMP, 1102 "Input array must be sparseness pattern (no values)"); 1103 } 1104 1105 if (power<=0) 1106 { 1107 SETERRQ( PETSC_ERR_SUP_SYS, 1108 "Power must be 1 or up"); 1109 } 1110 1111 // Copy input values 1112 retval.nrows = ncols; 1113 retval.ncols = nrows; 1114 retval.nnz = 0; 1115 retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS; 1116 retval.block_data = PETSC_FALSE; 1117 1118 // Allocate sparseness pattern 1119 ierr = spbas_allocate_pattern(&retval, in_matrix.values != NULL); 1120 CHKERRQ(ierr); 1121 1122 // Allocate marker array 1123 ierr = PetscMalloc(nrows * sizeof(PetscInt), &iwork); CHKERRQ(ierr); 1124 1125 // Erase the pattern for this row 1126 memset( (void *) iwork, 0, retval.nrows*sizeof(PetscInt)); 1127 1128 // Calculate marker values 1129 marker = 1; for (i=1; i<power; i++) {marker*=2;} 1130 1131 for (i=0; i<nrows; i++) 1132 { 1133 // Calculate the pattern for each row 1134 1135 nnz = in_matrix.row_nnz[i]; 1136 kend = i+in_matrix.icols[i][nnz-1]; 1137 if (maxmrk<=kend) {maxmrk=kend+1;} 1138 ierr = spbas_mark_row_power( iwork, i, &in_matrix, marker, 1139 i, maxmrk); CHKERRQ(ierr); 1140 1141 // Count the columns 1142 nnz = 0; 1143 for (j=i; j<maxmrk; j++) {nnz+= (iwork[j]!=0);} 1144 1145 // Allocate the column indices 1146 retval.row_nnz[i] = nnz; 1147 ierr = PetscMalloc(nnz*sizeof(PetscInt),&retval.icols[i]); CHKERRQ(ierr); 1148 1149 // Administrate the column indices 1150 inz = 0; 1151 for (j=i; j<maxmrk; j++) 1152 { 1153 if (iwork[j]) 1154 { 1155 retval.icols[i][inz] = j-i; 1156 inz++; 1157 iwork[j]=0; 1158 } 1159 } 1160 retval.nnz += nnz; 1161 }; 1162 1163 ierr = PetscFree(iwork); CHKERRQ(ierr); 1164 1165 *result = retval; 1166 PetscFunctionReturn(0); 1167 } 1168 1169 1170 1171 /* 1172 spbas_keep_upper: 1173 remove the lower part of the matrix: keep the upper part 1174 */ 1175 #undef __FUNCT__ 1176 #define __FUNCT__ "spbas_keep_upper" 1177 PetscErrorCode spbas_keep_upper( spbas_matrix * inout_matrix) 1178 { 1179 PetscInt i, j; 1180 PetscInt jstart; 1181 1182 PetscFunctionBegin; 1183 if (inout_matrix->block_data) 1184 { 1185 SETERRQ( PETSC_ERR_SUP_SYS, 1186 "Not yet for block data matrices\n"); 1187 } 1188 1189 for (i=0; i<inout_matrix->nrows; i++) 1190 { 1191 for (jstart=0; 1192 (jstart<inout_matrix->row_nnz[i]) && 1193 (inout_matrix->icols[i][jstart]<0); jstart++){} 1194 1195 if (jstart>0) 1196 { 1197 for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) 1198 { 1199 inout_matrix->icols[i][j] = 1200 inout_matrix->icols[i][j+jstart]; 1201 } 1202 1203 if (inout_matrix->values) 1204 { 1205 for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) 1206 { 1207 inout_matrix->values[i][j] = 1208 inout_matrix->values[i][j+jstart]; 1209 } 1210 } 1211 1212 inout_matrix->row_nnz[i] -= jstart; 1213 1214 inout_matrix->icols[i] = (PetscInt*) realloc( 1215 (void*) inout_matrix->icols[i], 1216 inout_matrix->row_nnz[i]*sizeof(PetscInt)); 1217 1218 if (inout_matrix->values) 1219 { 1220 inout_matrix->values[i] = (PetscScalar*) realloc( 1221 (void*) inout_matrix->values[i], 1222 inout_matrix->row_nnz[i]*sizeof(PetscScalar)); 1223 } 1224 inout_matrix->nnz -= jstart; 1225 } 1226 } 1227 PetscFunctionReturn(0); 1228 } 1229 1230