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", (PetscReal) nnz / (PetscReal) nrows); 422 423 424 // Clean up 425 ierr=PetscFree(isort); CHKERRQ(ierr); 426 ierr=PetscFree(used); CHKERRQ(ierr); 427 ierr=PetscFree(ipoint); CHKERRQ(ierr); 428 429 mem_compressed = spbas_memory_requirement( *B ); 430 *mem_reduction = 100.0 * (PetscScalar)(mem_orig-mem_compressed)/ 431 (PetscScalar) mem_orig; 432 PetscFunctionReturn(0); 433 } 434 435 /* 436 spbas_incomplete_cholesky 437 Incomplete Cholesky decomposition 438 */ 439 #include "spbas_cholesky.h" 440 441 /* 442 spbas_delete : de-allocate the arrays owned by this matrix 443 */ 444 #undef __FUNCT__ 445 #define __FUNCT__ "spbas_delete" 446 PetscErrorCode spbas_delete(spbas_matrix matrix) 447 { 448 PetscInt i; 449 PetscErrorCode ierr; 450 PetscFunctionBegin; 451 if (matrix.block_data) 452 { 453 ierr=PetscFree(matrix.alloc_icol); CHKERRQ(ierr) 454 if (matrix.values){ierr=PetscFree(matrix.alloc_val);CHKERRQ(ierr);} 455 } 456 else 457 { 458 for (i=0; i<matrix.nrows; i++) 459 { ierr=PetscFree(matrix.icols[i]); CHKERRQ(ierr);} 460 ierr = PetscFree(matrix.icols);CHKERRQ(ierr); 461 if (matrix.values) 462 { 463 for (i=0; i<matrix.nrows; i++) 464 { ierr=PetscFree(matrix.values[i]); CHKERRQ(ierr);} 465 } 466 } 467 468 ierr=PetscFree(matrix.row_nnz); CHKERRQ(ierr); 469 ierr=PetscFree(matrix.icols); CHKERRQ(ierr); 470 if (matrix.col_idx_type == SPBAS_OFFSET_ARRAY) 471 {ierr=PetscFree(matrix.icol0);CHKERRQ(ierr);} 472 if (matrix.values) 473 {ierr=PetscFree(matrix.values);CHKERRQ(ierr);} 474 475 PetscFunctionReturn(0); 476 } 477 478 /* 479 spbas_matrix_to_crs: 480 Convert an spbas_matrix to compessed row storage 481 */ 482 #undef __FUNCT__ 483 #define __FUNCT__ "spbas_matrix_to_crs" 484 PetscErrorCode spbas_matrix_to_crs(spbas_matrix matrix_A,MatScalar **val_out, PetscInt **irow_out, PetscInt **icol_out) 485 { 486 PetscInt nrows = matrix_A.nrows; 487 PetscInt nnz = matrix_A.nnz; 488 PetscInt i,j,r_nnz,i0; 489 PetscInt *irow; 490 PetscInt *icol; 491 PetscInt *icol_A; 492 MatScalar *val; 493 PetscScalar *val_A; 494 PetscInt col_idx_type = matrix_A.col_idx_type; 495 PetscTruth do_values = matrix_A.values ? PETSC_TRUE : PETSC_FALSE; 496 PetscErrorCode ierr; 497 498 PetscFunctionBegin; 499 ierr = PetscMalloc( sizeof(PetscInt) * (nrows+1), &irow); CHKERRQ(ierr); 500 ierr = PetscMalloc( sizeof(PetscInt) * nnz, &icol); CHKERRQ(ierr); 501 *icol_out = icol; *irow_out=irow; 502 if (do_values) 503 { 504 ierr = PetscMalloc( sizeof(MatScalar) * nnz, &val); CHKERRQ(ierr); 505 *val_out = val; *icol_out = icol; *irow_out=irow; 506 } 507 508 irow[0]=0; 509 for (i=0; i<nrows; i++) 510 { 511 r_nnz = matrix_A.row_nnz[i]; 512 i0 = irow[i]; 513 irow[i+1] = i0 + r_nnz; 514 icol_A = matrix_A.icols[i]; 515 516 if (do_values) 517 { 518 val_A = matrix_A.values[i]; 519 for (j=0; j<r_nnz; j++) 520 { 521 icol[i0+j] = icol_A[j]; 522 val[i0+j] = val_A[j]; 523 } 524 } 525 else 526 { 527 for (j=0; j<r_nnz; j++) { icol[i0+j] = icol_A[j]; } 528 } 529 530 if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 531 { 532 for (j=0; j<r_nnz; j++) { icol[i0+j] += i; } 533 } 534 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 535 { 536 i0 = matrix_A.icol0[i]; 537 for (j=0; j<r_nnz; j++) { icol[i0+j] += i0; } 538 } 539 540 } 541 542 PetscFunctionReturn(0); 543 } 544 545 /* 546 spbas_dump 547 Print the matrix in i,j,val-format 548 */ 549 #undef __FUNCT__ 550 #define __FUNCT__ "spbas_dump" 551 PetscErrorCode spbas_dump( const char *fname, spbas_matrix in_matrix) 552 { 553 FILE *outfile = fopen(fname, "w"); 554 PetscInt col_idx_type = in_matrix.col_idx_type; 555 PetscInt nrows=in_matrix.nrows; 556 PetscInt *icol; 557 PetscScalar *val; 558 PetscInt r_nnz, icol0; 559 PetscInt i,j; 560 PetscInt nwrite=1; 561 562 PetscFunctionBegin; 563 if (!outfile) 564 { 565 SETERRQ1( PETSC_ERR_FILE_OPEN, 566 "Error in spbas_dump: cannot open file '%s'\n",fname); 567 } 568 if (in_matrix.values) 569 { 570 for (i=0; i<nrows; i++) 571 { 572 icol = in_matrix.icols[i]; 573 val = in_matrix.values[i]; 574 r_nnz= in_matrix.row_nnz[i]; 575 if (col_idx_type == SPBAS_COLUMN_NUMBERS) 576 { 577 for (j=0; nwrite>0 && j<r_nnz; j++) 578 { nwrite = fprintf(outfile,"%d %d %e\n",i,icol[j],val[j]); } 579 } 580 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 581 { 582 for (j=0; nwrite>0 && j<r_nnz; j++) 583 { nwrite = fprintf(outfile,"%d %d %e\n",i,i+icol[j],val[j]);} 584 } 585 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 586 { 587 icol0 = in_matrix.icol0[i]; 588 for (j=0; nwrite>0 && j<r_nnz; j++) 589 { nwrite = fprintf(outfile,"%d %d %e\n",i,icol0+icol[j],val[j]); } 590 } 591 592 } 593 } 594 else 595 { 596 for (i=0; i<nrows; i++) 597 { 598 icol = in_matrix.icols[i]; 599 r_nnz= in_matrix.row_nnz[i]; 600 nwrite = 2; 601 if (col_idx_type == SPBAS_COLUMN_NUMBERS) 602 { 603 for (j=0; nwrite>0 && j<r_nnz; j++) 604 { nwrite = fprintf(outfile,"%d %d\n",i,icol[j]); } 605 } 606 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 607 { 608 for (j=0; nwrite>0 && j<r_nnz; j++) 609 { nwrite = fprintf(outfile,"%d %d\n",i,i+icol[j]); } 610 } 611 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 612 { 613 icol0 = in_matrix.icol0[i]; 614 for (j=0; nwrite>0 && j<r_nnz; j++) 615 { nwrite = fprintf(outfile,"%d %d\n",i,icol0+icol[j]); } 616 } 617 } 618 } 619 if (nwrite<0) 620 { 621 SETERRQ1(PETSC_ERR_FILE_WRITE, 622 "Error in spbas_dump: cannot write to file '%s'\n",fname); 623 } 624 fclose(outfile); 625 PetscFunctionReturn(0); 626 } 627 628 /* 629 spbas_transpose 630 return the transpose of a matrix 631 */ 632 #undef __FUNCT__ 633 #define __FUNCT__ "spbas_transpose" 634 PetscErrorCode spbas_transpose( spbas_matrix in_matrix, spbas_matrix * result) 635 { 636 PetscInt col_idx_type = in_matrix.col_idx_type; 637 PetscInt nnz = in_matrix.nnz; 638 PetscInt ncols = in_matrix.nrows; 639 PetscInt nrows = in_matrix.ncols; 640 PetscInt i,j,k; 641 PetscInt r_nnz; 642 PetscInt *irow; 643 PetscInt icol0; 644 PetscScalar * val; 645 PetscErrorCode ierr; 646 647 PetscFunctionBegin; 648 649 // Copy input values 650 result->nrows = nrows; 651 result->ncols = ncols; 652 result->nnz = nnz; 653 result->col_idx_type = SPBAS_COLUMN_NUMBERS; 654 result->block_data = PETSC_TRUE; 655 656 // Allocate sparseness pattern 657 ierr = spbas_allocate_pattern(result, in_matrix.values ? PETSC_TRUE : PETSC_FALSE); 658 CHKERRQ(ierr); 659 660 // Count the number of nonzeros in each row 661 for (i = 0; i<nrows; i++) { result->row_nnz[i] = 0; } 662 663 for (i=0; i<ncols; i++) 664 { 665 r_nnz = in_matrix.row_nnz[i]; 666 irow = in_matrix.icols[i]; 667 if (col_idx_type == SPBAS_COLUMN_NUMBERS) 668 { 669 for (j=0; j<r_nnz; j++) { result->row_nnz[irow[j]]++; } 670 } 671 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) 672 { 673 for (j=0; j<r_nnz; j++) { result->row_nnz[i+irow[j]]++; } 674 } 675 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 676 { 677 icol0=in_matrix.icol0[i]; 678 for (j=0; j<r_nnz; j++) { result->row_nnz[icol0+irow[j]]++; } 679 } 680 } 681 682 // Set the pointers to the data 683 ierr = spbas_allocate_data(result); CHKERRQ(ierr); 684 685 // Reset the number of nonzeros in each row 686 for (i = 0; i<nrows; i++) { result->row_nnz[i] = 0; } 687 688 // Fill the data arrays 689 if (in_matrix.values) 690 { 691 for (i=0; i<ncols; i++) 692 { 693 r_nnz = in_matrix.row_nnz[i]; 694 irow = in_matrix.icols[i]; 695 val = in_matrix.values[i]; 696 697 if (col_idx_type == SPBAS_COLUMN_NUMBERS) {icol0=0;} 698 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {icol0=i;} 699 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 700 {icol0=in_matrix.icol0[i];} 701 for (j=0; j<r_nnz; j++) 702 { 703 k = icol0 + irow[j]; 704 result->icols[k][result->row_nnz[k]] = i; 705 result->values[k][result->row_nnz[k]] = val[j]; 706 result->row_nnz[k]++; 707 } 708 } 709 } 710 else 711 { 712 for (i=0; i<ncols; i++) 713 { 714 r_nnz = in_matrix.row_nnz[i]; 715 irow = in_matrix.icols[i]; 716 717 if (col_idx_type == SPBAS_COLUMN_NUMBERS) {icol0=0;} 718 else if (col_idx_type == SPBAS_DIAGONAL_OFFSETS) {icol0=i;} 719 else if (col_idx_type == SPBAS_OFFSET_ARRAY) 720 {icol0=in_matrix.icol0[i];} 721 722 for (j=0; j<r_nnz; j++) 723 { 724 k = icol0 + irow[j]; 725 result->icols[k][result->row_nnz[k]] = i; 726 result->row_nnz[k]++; 727 } 728 } 729 } 730 731 // Set return value 732 PetscFunctionReturn(0); 733 } 734 735 /* 736 spbas_mergesort 737 738 mergesort for an array of intergers and an array of associated 739 reals 740 741 on output, icol[0..nnz-1] is increasing; 742 val[0..nnz-1] has undergone the same permutation as icol 743 744 NB: val may be NULL: in that case, only the integers are sorted 745 746 */ 747 #undef __FUNCT__ 748 #define __FUNCT__ "spbas_mergesort" 749 PetscErrorCode spbas_mergesort(PetscInt nnz, PetscInt *icol, PetscScalar *val) 750 { 751 PetscInt istep; // Chunk-sizes of already sorted parts of arrays 752 PetscInt i, i1, i2; // Loop counters for (partly) sorted arrays 753 PetscInt istart, i1end, i2end; // start of newly sorted array part, end of both parts 754 PetscInt *ialloc; // Allocated arrays 755 PetscScalar *valloc=NULL; 756 PetscInt *iswap; // auxiliary pointers for swapping 757 PetscScalar *vswap; 758 PetscInt *ihlp1; // Pointers to new version of arrays, 759 PetscScalar *vhlp1=NULL; // (arrays under construction) 760 PetscInt *ihlp2; // Pointers to previous version of arrays, 761 PetscScalar *vhlp2=NULL; 762 PetscErrorCode ierr; 763 764 // Create arrays 765 ierr = PetscMalloc(nnz*sizeof(PetscInt),&ialloc); CHKERRQ(ierr); 766 ihlp1 = ialloc; 767 ihlp2 = icol; 768 769 if (val) 770 { 771 ierr = PetscMalloc(nnz*sizeof(PetscScalar),&valloc); CHKERRQ(ierr); 772 vhlp1 = valloc; 773 vhlp2 = val; 774 } 775 776 777 // Sorted array chunks are first 1 long, and increase until they are the 778 // complete array 779 for (istep=1; istep<nnz; istep*=2) 780 { 781 // 782 // Combine sorted parts 783 // istart:istart+istep-1 and istart+istep-1:istart+2*istep-1 784 // of ihlp2 and vhlp2 785 // 786 // into one sorted part 787 // istart:istart+2*istep-1 788 // of ihlp1 and vhlp1 789 // 790 for (istart=0; istart<nnz; istart+=2*istep) 791 { 792 793 // Set counters and bound array part endings 794 i1=istart; i1end = i1+istep; if (i1end>nnz) {i1end=nnz;} 795 i2=istart+istep; i2end = i2+istep; if (i2end>nnz) {i2end=nnz;} 796 797 // Merge the two array parts 798 if (val) 799 { 800 for (i=istart; i<i2end; i++) 801 { 802 if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) 803 { 804 ihlp1[i] = ihlp2[i1]; 805 vhlp1[i] = vhlp2[i1]; 806 i1++; 807 } 808 else if (i2<i2end ) 809 { 810 ihlp1[i] = ihlp2[i2]; 811 vhlp1[i] = vhlp2[i2]; 812 i2++; 813 } 814 else 815 { 816 ihlp1[i] = ihlp2[i1]; 817 vhlp1[i] = vhlp2[i1]; 818 i1++; 819 } 820 } 821 } 822 else 823 { 824 for (i=istart; i<i2end; i++) 825 { 826 if (i1<i1end && i2<i2end && ihlp2[i1] < ihlp2[i2]) 827 { 828 ihlp1[i] = ihlp2[i1]; 829 i1++; 830 } 831 else if (i2<i2end ) 832 { 833 ihlp1[i] = ihlp2[i2]; 834 i2++; 835 } 836 else 837 { 838 ihlp1[i] = ihlp2[i1]; 839 i1++; 840 } 841 } 842 } 843 } 844 845 // Swap the two array sets 846 iswap = ihlp2; ihlp2 = ihlp1; ihlp1 = iswap; 847 vswap = vhlp2; vhlp2 = vhlp1; vhlp1 = vswap; 848 } 849 850 // Copy one more time in case the sorted arrays are the temporary ones 851 if (ihlp2 != icol) 852 { 853 for (i=0; i<nnz; i++) { icol[i] = ihlp2[i]; } 854 if (val) 855 { 856 for (i=0; i<nnz; i++) { val[i] = vhlp2[i]; } 857 } 858 } 859 860 // Remove help arrays 861 ierr = PetscFree(ialloc); CHKERRQ(ierr); 862 if(val){ierr = PetscFree(valloc); CHKERRQ(ierr);} 863 PetscFunctionReturn(0); 864 } 865 866 /* 867 spbas_apply_reordering_rows: 868 apply the given reordering to the rows: matrix_A = matrix_A(perm,:); 869 */ 870 #undef __FUNCT__ 871 #define __FUNCT__ "spbas_apply_reordering_rows" 872 PetscErrorCode spbas_apply_reordering_rows(spbas_matrix *matrix_A, const PetscInt *permutation) 873 { 874 PetscInt i,j,ip; 875 PetscInt nrows=matrix_A->nrows; 876 PetscInt * row_nnz; 877 PetscInt **icols; 878 PetscTruth do_values = matrix_A->values ? PETSC_TRUE : PETSC_FALSE; 879 PetscScalar **vals=NULL; 880 PetscErrorCode ierr; 881 882 PetscFunctionBegin; 883 if (matrix_A->col_idx_type != SPBAS_DIAGONAL_OFFSETS) 884 { 885 SETERRQ( PETSC_ERR_SUP_SYS, 886 "must have diagonal offsets in pattern\n"); 887 } 888 889 if (do_values) 890 { 891 ierr = PetscMalloc( sizeof(PetscScalar*)*nrows, &vals); CHKERRQ(ierr); 892 } 893 ierr = PetscMalloc( sizeof(PetscInt)*nrows, &row_nnz); CHKERRQ(ierr); 894 ierr = PetscMalloc( sizeof(PetscInt*)*nrows, &icols); CHKERRQ(ierr); 895 896 for (i=0; i<nrows;i++) 897 { 898 ip = permutation[i]; 899 if (do_values) {vals[i] = matrix_A->values[ip];} 900 icols[i] = matrix_A->icols[ip]; 901 row_nnz[i] = matrix_A->row_nnz[ip]; 902 for (j=0; j<row_nnz[i]; j++) { icols[i][j] += ip-i; } 903 } 904 905 if (do_values){ ierr = PetscFree(matrix_A->values); CHKERRQ(ierr);} 906 ierr = PetscFree(matrix_A->icols); CHKERRQ(ierr); 907 ierr = PetscFree(matrix_A->row_nnz); CHKERRQ(ierr); 908 909 if (do_values) { matrix_A->values = vals; } 910 matrix_A->icols = icols; 911 matrix_A->row_nnz = row_nnz; 912 913 PetscFunctionReturn(0); 914 } 915 916 917 /* 918 spbas_apply_reordering_cols: 919 apply the given reordering to the columns: matrix_A(:,perm) = matrix_A; 920 */ 921 #undef __FUNCT__ 922 #define __FUNCT__ "spbas_apply_reordering_cols" 923 PetscErrorCode spbas_apply_reordering_cols( spbas_matrix *matrix_A,const PetscInt *permutation) 924 { 925 PetscInt i,j; 926 PetscInt nrows=matrix_A->nrows; 927 PetscInt row_nnz; 928 PetscInt *icols; 929 PetscTruth do_values = matrix_A->values ? PETSC_TRUE : PETSC_FALSE; 930 PetscScalar *vals=NULL; 931 932 PetscErrorCode ierr; 933 934 PetscFunctionBegin; 935 936 if (matrix_A->col_idx_type != SPBAS_DIAGONAL_OFFSETS) 937 { 938 SETERRQ( PETSC_ERR_SUP_SYS, 939 "must have diagonal offsets in pattern\n"); 940 } 941 942 for (i=0; i<nrows;i++) 943 { 944 icols = matrix_A->icols[i]; 945 row_nnz = matrix_A->row_nnz[i]; 946 if (do_values) { vals = matrix_A->values[i]; } 947 948 for (j=0; j<row_nnz; j++) 949 { 950 icols[j] = permutation[i+icols[j]]-i; 951 } 952 ierr = spbas_mergesort(row_nnz, icols, vals); CHKERRQ(ierr); 953 } 954 955 PetscFunctionReturn(0); 956 } 957 958 /* 959 spbas_apply_reordering: 960 apply the given reordering: matrix_A(perm,perm) = matrix_A; 961 */ 962 #undef __FUNCT__ 963 #define __FUNCT__ "spbas_apply_reordering" 964 PetscErrorCode spbas_apply_reordering( spbas_matrix *matrix_A, const PetscInt *permutation, const PetscInt * inv_perm) 965 { 966 PetscErrorCode ierr; 967 PetscFunctionBegin; 968 ierr = spbas_apply_reordering_rows( matrix_A, inv_perm); CHKERRQ(ierr); 969 ierr = spbas_apply_reordering_cols( matrix_A, permutation); CHKERRQ(ierr); 970 PetscFunctionReturn(0); 971 } 972 973 #undef __FUNCT__ 974 #define __FUNCT__ "spbas_pattern_only" 975 PetscErrorCode spbas_pattern_only(PetscInt nrows, PetscInt ncols, PetscInt *ai, PetscInt *aj, spbas_matrix * result) 976 { 977 spbas_matrix retval; 978 PetscInt i, j, i0, r_nnz; 979 PetscErrorCode ierr; 980 981 PetscFunctionBegin; 982 983 // Copy input values 984 retval.nrows = nrows; 985 retval.ncols = ncols; 986 retval.nnz = ai[nrows]; 987 988 retval.block_data = PETSC_TRUE; 989 retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS; 990 991 // Allocate output matrix 992 ierr = spbas_allocate_pattern(&retval, PETSC_FALSE); CHKERRQ(ierr); 993 for (i=0; i<nrows; i++) {retval.row_nnz[i] = ai[i+1]-ai[i];} 994 ierr = spbas_allocate_data(&retval); CHKERRQ(ierr); 995 // Copy the structure 996 for (i = 0; i<retval.nrows; i++) 997 { 998 i0 = ai[i]; 999 r_nnz = ai[i+1]-i0; 1000 1001 for (j=0; j<r_nnz; j++) 1002 { 1003 retval.icols[i][j] = aj[i0+j]-i; 1004 } 1005 } 1006 1007 // Set return value 1008 *result = retval; 1009 PetscFunctionReturn(0); 1010 } 1011 1012 1013 /* 1014 spbas_mark_row_power: 1015 Mark the columns in row 'row' which are nonzero in 1016 matrix^2log(marker). 1017 */ 1018 #undef __FUNCT__ 1019 #define __FUNCT__ "spbas_mark_row_power" 1020 PetscErrorCode spbas_mark_row_power( 1021 PetscInt *iwork, // marker-vector 1022 PetscInt row, // row for which the columns are marked 1023 spbas_matrix * in_matrix, // matrix for which the power is being 1024 // calculated 1025 PetscInt marker, // marker-value: 2^power 1026 PetscInt minmrk, // lower bound for marked points 1027 PetscInt maxmrk) // upper bound for marked points 1028 { 1029 PetscErrorCode ierr; 1030 PetscInt i,j, nnz; 1031 1032 PetscFunctionBegin; 1033 nnz = in_matrix->row_nnz[row]; 1034 1035 // For higher powers, call this function recursively 1036 if (marker>1) 1037 { 1038 for (i=0; i<nnz;i++) 1039 { 1040 j = row + in_matrix->icols[row][i]; 1041 if (minmrk<=j && j<maxmrk && iwork[j] < marker ) 1042 { 1043 ierr = spbas_mark_row_power( iwork, row + in_matrix->icols[row][i], 1044 in_matrix, marker/2,minmrk,maxmrk); 1045 CHKERRQ(ierr); 1046 iwork[j] |= marker; 1047 } 1048 } 1049 } 1050 else 1051 { 1052 // Mark the columns reached 1053 for (i=0; i<nnz;i++) 1054 { 1055 j = row + in_matrix->icols[row][i]; 1056 if (minmrk<=j && j<maxmrk) { iwork[j] |= 1; } 1057 } 1058 } 1059 1060 PetscFunctionReturn(0); 1061 } 1062 1063 1064 /* 1065 spbas_power 1066 Calculate sparseness patterns for incomplete Cholesky decompositions 1067 of a given order: (almost) all nonzeros of the matrix^(order+1) which 1068 are inside the band width are found and stored in the output sparseness 1069 pattern. 1070 */ 1071 #undef __FUNCT__ 1072 #define __FUNCT__ "spbas_power" 1073 PetscErrorCode spbas_power (spbas_matrix in_matrix,PetscInt power, spbas_matrix * result) 1074 { 1075 spbas_matrix retval; 1076 PetscInt nrows = in_matrix.nrows; 1077 PetscInt ncols = in_matrix.ncols; 1078 PetscInt i, j, kend; 1079 PetscInt nnz, inz; 1080 PetscInt *iwork; 1081 PetscInt marker; 1082 PetscInt maxmrk=0; 1083 PetscErrorCode ierr; 1084 1085 PetscFunctionBegin; 1086 if (in_matrix.col_idx_type != SPBAS_DIAGONAL_OFFSETS) 1087 { 1088 SETERRQ( PETSC_ERR_SUP_SYS, 1089 "must have diagonal offsets in pattern\n"); 1090 } 1091 1092 if (ncols != nrows) 1093 { 1094 SETERRQ( PETSC_ERR_ARG_INCOMP, 1095 "Dimension error\n"); 1096 } 1097 1098 if (in_matrix.values) 1099 { 1100 SETERRQ( PETSC_ERR_ARG_INCOMP, 1101 "Input array must be sparseness pattern (no values)"); 1102 } 1103 1104 if (power<=0) 1105 { 1106 SETERRQ( PETSC_ERR_SUP_SYS, 1107 "Power must be 1 or up"); 1108 } 1109 1110 // Copy input values 1111 retval.nrows = ncols; 1112 retval.ncols = nrows; 1113 retval.nnz = 0; 1114 retval.col_idx_type = SPBAS_DIAGONAL_OFFSETS; 1115 retval.block_data = PETSC_FALSE; 1116 1117 // Allocate sparseness pattern 1118 ierr = spbas_allocate_pattern(&retval, in_matrix.values ? PETSC_TRUE : PETSC_FALSE); 1119 CHKERRQ(ierr); 1120 1121 // Allocate marker array 1122 ierr = PetscMalloc(nrows * sizeof(PetscInt), &iwork); CHKERRQ(ierr); 1123 1124 // Erase the pattern for this row 1125 memset( (void *) iwork, 0, retval.nrows*sizeof(PetscInt)); 1126 1127 // Calculate marker values 1128 marker = 1; for (i=1; i<power; i++) {marker*=2;} 1129 1130 for (i=0; i<nrows; i++) 1131 { 1132 // Calculate the pattern for each row 1133 1134 nnz = in_matrix.row_nnz[i]; 1135 kend = i+in_matrix.icols[i][nnz-1]; 1136 if (maxmrk<=kend) {maxmrk=kend+1;} 1137 ierr = spbas_mark_row_power( iwork, i, &in_matrix, marker, 1138 i, maxmrk); CHKERRQ(ierr); 1139 1140 // Count the columns 1141 nnz = 0; 1142 for (j=i; j<maxmrk; j++) {nnz+= (iwork[j]!=0);} 1143 1144 // Allocate the column indices 1145 retval.row_nnz[i] = nnz; 1146 ierr = PetscMalloc(nnz*sizeof(PetscInt),&retval.icols[i]); CHKERRQ(ierr); 1147 1148 // Administrate the column indices 1149 inz = 0; 1150 for (j=i; j<maxmrk; j++) 1151 { 1152 if (iwork[j]) 1153 { 1154 retval.icols[i][inz] = j-i; 1155 inz++; 1156 iwork[j]=0; 1157 } 1158 } 1159 retval.nnz += nnz; 1160 }; 1161 1162 ierr = PetscFree(iwork); CHKERRQ(ierr); 1163 1164 *result = retval; 1165 PetscFunctionReturn(0); 1166 } 1167 1168 1169 1170 /* 1171 spbas_keep_upper: 1172 remove the lower part of the matrix: keep the upper part 1173 */ 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "spbas_keep_upper" 1176 PetscErrorCode spbas_keep_upper( spbas_matrix * inout_matrix) 1177 { 1178 PetscInt i, j; 1179 PetscInt jstart; 1180 1181 PetscFunctionBegin; 1182 if (inout_matrix->block_data) 1183 { 1184 SETERRQ( PETSC_ERR_SUP_SYS, 1185 "Not yet for block data matrices\n"); 1186 } 1187 1188 for (i=0; i<inout_matrix->nrows; i++) 1189 { 1190 for (jstart=0; 1191 (jstart<inout_matrix->row_nnz[i]) && 1192 (inout_matrix->icols[i][jstart]<0); jstart++){} 1193 1194 if (jstart>0) 1195 { 1196 for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) 1197 { 1198 inout_matrix->icols[i][j] = 1199 inout_matrix->icols[i][j+jstart]; 1200 } 1201 1202 if (inout_matrix->values) 1203 { 1204 for (j=0; j<inout_matrix->row_nnz[i]-jstart; j++) 1205 { 1206 inout_matrix->values[i][j] = 1207 inout_matrix->values[i][j+jstart]; 1208 } 1209 } 1210 1211 inout_matrix->row_nnz[i] -= jstart; 1212 1213 inout_matrix->icols[i] = (PetscInt*) realloc( 1214 (void*) inout_matrix->icols[i], 1215 inout_matrix->row_nnz[i]*sizeof(PetscInt)); 1216 1217 if (inout_matrix->values) 1218 { 1219 inout_matrix->values[i] = (PetscScalar*) realloc( 1220 (void*) inout_matrix->values[i], 1221 inout_matrix->row_nnz[i]*sizeof(PetscScalar)); 1222 } 1223 inout_matrix->nnz -= jstart; 1224 } 1225 } 1226 PetscFunctionReturn(0); 1227 } 1228 1229