1 /*$Id: baij.c,v 1.245 2001/08/07 03:02:55 balay Exp bsmith $*/ 2 3 /* 4 Defines the basic matrix operations for the BAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "src/mat/impls/baij/seq/baij.h" 8 #include "src/vec/vecimpl.h" 9 #include "src/inline/spops.h" 10 #include "petscsys.h" /*I "petscmat.h" I*/ 11 12 /* 13 Special version for Fun3d sequential benchmark 14 */ 15 #if defined(PETSC_HAVE_FORTRAN_CAPS) 16 #define matsetvaluesblocked4_ MATSETVALUESBLOCKED4 17 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 18 #define matsetvaluesblocked4_ matsetvaluesblocked4 19 #endif 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "matsetvaluesblocked4_" 23 void matsetvaluesblocked4_(Mat *AA,int *mm,int *im,int *nn,int *in,MatScalar *v,InsertMode *is,int *err) 24 { 25 Mat A = *AA; 26 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 27 int *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,m = *mm,n = *nn; 28 int *imax=a->imax,*ai=a->i,*ailen=a->ilen; 29 int *aj=a->j,nonew=a->nonew,stepval,ierr; 30 MatScalar *value = v,*ap,*aa = a->a,*bap; 31 32 PetscFunctionBegin; 33 stepval = (n-1)*4; 34 for (k=0; k<m; k++) { /* loop over added rows */ 35 row = im[k]; 36 rp = aj + ai[row]; 37 ap = aa + 16*ai[row]; 38 rmax = imax[row]; 39 nrow = ailen[row]; 40 low = 0; 41 for (l=0; l<n; l++) { /* loop over added columns */ 42 col = in[l]; 43 value = v + k*(stepval+4)*4 + l*4; 44 low = 0; high = nrow; 45 while (high-low > 7) { 46 t = (low+high)/2; 47 if (rp[t] > col) high = t; 48 else low = t; 49 } 50 for (i=low; i<high; i++) { 51 if (rp[i] > col) break; 52 if (rp[i] == col) { 53 bap = ap + 16*i; 54 for (ii=0; ii<4; ii++,value+=stepval) { 55 for (jj=ii; jj<16; jj+=4) { 56 bap[jj] += *value++; 57 } 58 } 59 goto noinsert2; 60 } 61 } 62 N = nrow++ - 1; 63 /* shift up all the later entries in this row */ 64 for (ii=N; ii>=i; ii--) { 65 rp[ii+1] = rp[ii]; 66 ierr = PetscMemcpy(ap+16*(ii+1),ap+16*(ii),16*sizeof(MatScalar)); 67 } 68 if (N >= i) { 69 ierr = PetscMemzero(ap+16*i,16*sizeof(MatScalar)); 70 } 71 rp[i] = col; 72 bap = ap + 16*i; 73 for (ii=0; ii<4; ii++,value+=stepval) { 74 for (jj=ii; jj<16; jj+=4) { 75 bap[jj] = *value++; 76 } 77 } 78 noinsert2:; 79 low = i; 80 } 81 ailen[row] = nrow; 82 } 83 } 84 85 #if defined(PETSC_HAVE_FORTRAN_CAPS) 86 #define matsetvalues4_ MATSETVALUES4 87 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 88 #define matsetvalues4_ matsetvalues4 89 #endif 90 91 #undef __FUNCT__ 92 #define __FUNCT__ "MatSetValues4_" 93 int matsetvalues4_(Mat *AA,int *mm,int *im,int *nn,int *in,PetscScalar *v,InsertMode *is,int *err) 94 { 95 Mat A = *AA; 96 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 97 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted,n = *nn,m = *mm; 98 int *imax=a->imax,*ai=a->i,*ailen=a->ilen; 99 int *aj=a->j,brow,bcol; 100 int ridx,cidx,ierr; 101 MatScalar *ap,value,*aa=a->a,*bap; 102 103 PetscFunctionBegin; 104 for (k=0; k<m; k++) { /* loop over added rows */ 105 row = im[k]; brow = row/4; 106 rp = aj + ai[brow]; 107 ap = aa + 16*ai[brow]; 108 rmax = imax[brow]; 109 nrow = ailen[brow]; 110 low = 0; 111 for (l=0; l<n; l++) { /* loop over added columns */ 112 col = in[l]; bcol = col/4; 113 ridx = row % 4; cidx = col % 4; 114 value = v[l + k*n]; 115 low = 0; high = nrow; 116 while (high-low > 7) { 117 t = (low+high)/2; 118 if (rp[t] > bcol) high = t; 119 else low = t; 120 } 121 for (i=low; i<high; i++) { 122 if (rp[i] > bcol) break; 123 if (rp[i] == bcol) { 124 bap = ap + 16*i + 4*cidx + ridx; 125 *bap += value; 126 goto noinsert1; 127 } 128 } 129 N = nrow++ - 1; 130 /* shift up all the later entries in this row */ 131 for (ii=N; ii>=i; ii--) { 132 rp[ii+1] = rp[ii]; 133 ierr = PetscMemcpy(ap+16*(ii+1),ap+16*(ii),16*sizeof(MatScalar));CHKERRQ(ierr); 134 } 135 if (N>=i) { 136 ierr = PetscMemzero(ap+16*i,16*sizeof(MatScalar));CHKERRQ(ierr); 137 } 138 rp[i] = bcol; 139 ap[16*i + 4*cidx + ridx] = value; 140 noinsert1:; 141 low = i; 142 } 143 ailen[brow] = nrow; 144 } 145 } 146 147 /* UGLY, ugly, ugly 148 When MatScalar == PetscScalar the function MatSetValuesBlocked_SeqBAIJ_MatScalar() does 149 not exist. Otherwise ..._MatScalar() takes matrix dlements in single precision and 150 inserts them into the single precision data structure. The function MatSetValuesBlocked_SeqBAIJ() 151 converts the entries into single precision and then calls ..._MatScalar() to put them 152 into the single precision data structures. 153 */ 154 #if defined(PETSC_USE_MAT_SINGLE) 155 EXTERN int MatSetValuesBlocked_SeqBAIJ_MatScalar(Mat,int,int*,int,int*,MatScalar*,InsertMode); 156 #else 157 #define MatSetValuesBlocked_SeqBAIJ_MatScalar MatSetValuesBlocked_SeqBAIJ 158 #endif 159 #if defined(PETSC_HAVE_DSCPACK) 160 EXTERN int MatUseDSCPACK_MPIBAIJ(Mat); 161 #endif 162 163 #define CHUNKSIZE 10 164 165 /* 166 Checks for missing diagonals 167 */ 168 #undef __FUNCT__ 169 #define __FUNCT__ "MatMissingDiagonal_SeqBAIJ" 170 int MatMissingDiagonal_SeqBAIJ(Mat A) 171 { 172 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 173 int *diag,*jj = a->j,i,ierr; 174 175 PetscFunctionBegin; 176 ierr = MatMarkDiagonal_SeqBAIJ(A);CHKERRQ(ierr); 177 diag = a->diag; 178 for (i=0; i<a->mbs; i++) { 179 if (jj[diag[i]] != i) { 180 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 181 } 182 } 183 PetscFunctionReturn(0); 184 } 185 186 #undef __FUNCT__ 187 #define __FUNCT__ "MatMarkDiagonal_SeqBAIJ" 188 int MatMarkDiagonal_SeqBAIJ(Mat A) 189 { 190 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 191 int i,j,*diag,m = a->mbs,ierr; 192 193 PetscFunctionBegin; 194 if (a->diag) PetscFunctionReturn(0); 195 196 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 197 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 198 for (i=0; i<m; i++) { 199 diag[i] = a->i[i+1]; 200 for (j=a->i[i]; j<a->i[i+1]; j++) { 201 if (a->j[j] == i) { 202 diag[i] = j; 203 break; 204 } 205 } 206 } 207 a->diag = diag; 208 PetscFunctionReturn(0); 209 } 210 211 212 EXTERN int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "MatGetRowIJ_SeqBAIJ" 216 static int MatGetRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 217 { 218 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 219 int ierr,n = a->mbs,i; 220 221 PetscFunctionBegin; 222 *nn = n; 223 if (!ia) PetscFunctionReturn(0); 224 if (symmetric) { 225 ierr = MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 226 } else if (oshift == 1) { 227 /* temporarily add 1 to i and j indices */ 228 int nz = a->i[n]; 229 for (i=0; i<nz; i++) a->j[i]++; 230 for (i=0; i<n+1; i++) a->i[i]++; 231 *ia = a->i; *ja = a->j; 232 } else { 233 *ia = a->i; *ja = a->j; 234 } 235 236 PetscFunctionReturn(0); 237 } 238 239 #undef __FUNCT__ 240 #define __FUNCT__ "MatRestoreRowIJ_SeqBAIJ" 241 static int MatRestoreRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 242 { 243 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 244 int i,n = a->mbs,ierr; 245 246 PetscFunctionBegin; 247 if (!ia) PetscFunctionReturn(0); 248 if (symmetric) { 249 ierr = PetscFree(*ia);CHKERRQ(ierr); 250 ierr = PetscFree(*ja);CHKERRQ(ierr); 251 } else if (oshift == 1) { 252 int nz = a->i[n]-1; 253 for (i=0; i<nz; i++) a->j[i]--; 254 for (i=0; i<n+1; i++) a->i[i]--; 255 } 256 PetscFunctionReturn(0); 257 } 258 259 #undef __FUNCT__ 260 #define __FUNCT__ "MatGetBlockSize_SeqBAIJ" 261 int MatGetBlockSize_SeqBAIJ(Mat mat,int *bs) 262 { 263 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ*)mat->data; 264 265 PetscFunctionBegin; 266 *bs = baij->bs; 267 PetscFunctionReturn(0); 268 } 269 270 271 #undef __FUNCT__ 272 #define __FUNCT__ "MatDestroy_SeqBAIJ" 273 int MatDestroy_SeqBAIJ(Mat A) 274 { 275 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 276 int ierr; 277 278 PetscFunctionBegin; 279 #if defined(PETSC_USE_LOG) 280 PetscLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",A->m,A->n,a->nz); 281 #endif 282 ierr = PetscFree(a->a);CHKERRQ(ierr); 283 if (!a->singlemalloc) { 284 ierr = PetscFree(a->i);CHKERRQ(ierr); 285 ierr = PetscFree(a->j);CHKERRQ(ierr); 286 } 287 if (a->row) { 288 ierr = ISDestroy(a->row);CHKERRQ(ierr); 289 } 290 if (a->col) { 291 ierr = ISDestroy(a->col);CHKERRQ(ierr); 292 } 293 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 294 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 295 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 296 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 297 if (a->mult_work) {ierr = PetscFree(a->mult_work);CHKERRQ(ierr);} 298 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 299 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 300 #if defined(PETSC_USE_MAT_SINGLE) 301 if (a->setvaluescopy) {ierr = PetscFree(a->setvaluescopy);CHKERRQ(ierr);} 302 #endif 303 ierr = PetscFree(a);CHKERRQ(ierr); 304 PetscFunctionReturn(0); 305 } 306 307 #undef __FUNCT__ 308 #define __FUNCT__ "MatSetOption_SeqBAIJ" 309 int MatSetOption_SeqBAIJ(Mat A,MatOption op) 310 { 311 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 312 313 PetscFunctionBegin; 314 switch (op) { 315 case MAT_ROW_ORIENTED: 316 a->roworiented = PETSC_TRUE; 317 break; 318 case MAT_COLUMN_ORIENTED: 319 a->roworiented = PETSC_FALSE; 320 break; 321 case MAT_COLUMNS_SORTED: 322 a->sorted = PETSC_TRUE; 323 break; 324 case MAT_COLUMNS_UNSORTED: 325 a->sorted = PETSC_FALSE; 326 break; 327 case MAT_KEEP_ZEROED_ROWS: 328 a->keepzeroedrows = PETSC_TRUE; 329 break; 330 case MAT_NO_NEW_NONZERO_LOCATIONS: 331 a->nonew = 1; 332 break; 333 case MAT_NEW_NONZERO_LOCATION_ERR: 334 a->nonew = -1; 335 break; 336 case MAT_NEW_NONZERO_ALLOCATION_ERR: 337 a->nonew = -2; 338 break; 339 case MAT_YES_NEW_NONZERO_LOCATIONS: 340 a->nonew = 0; 341 break; 342 case MAT_ROWS_SORTED: 343 case MAT_ROWS_UNSORTED: 344 case MAT_YES_NEW_DIAGONALS: 345 case MAT_IGNORE_OFF_PROC_ENTRIES: 346 case MAT_USE_HASH_TABLE: 347 PetscLogInfo(A,"MatSetOption_SeqBAIJ:Option ignored\n"); 348 break; 349 case MAT_NO_NEW_DIAGONALS: 350 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 351 case MAT_USE_SINGLE_PRECISION_SOLVES: 352 if (a->bs==4) { 353 PetscTruth sse_enabled_local,sse_enabled_global; 354 int ierr; 355 356 sse_enabled_local = PETSC_FALSE; 357 sse_enabled_global = PETSC_FALSE; 358 359 ierr = PetscSSEIsEnabled(A->comm,&sse_enabled_local,&sse_enabled_global);CHKERRQ(ierr); 360 #if defined(PETSC_HAVE_SSE) 361 if (sse_enabled_local) { 362 a->single_precision_solves = PETSC_TRUE; 363 A->ops->solve = MatSolve_SeqBAIJ_Update; 364 A->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_Update; 365 PetscLogInfo(A,"MatSetOption_SeqBAIJ:Option MAT_USE_SINGLE_PRECISION_SOLVES set\n"); 366 break; 367 } else { 368 PetscLogInfo(A,"MatSetOption_SeqBAIJ:Option MAT_USE_SINGLE_PRECISION_SOLVES ignored\n"); 369 } 370 #else 371 PetscLogInfo(A,"MatSetOption_SeqBAIJ:Option MAT_USE_SINGLE_PRECISION_SOLVES ignored\n"); 372 #endif 373 } 374 break; 375 default: 376 SETERRQ(PETSC_ERR_SUP,"unknown option"); 377 } 378 PetscFunctionReturn(0); 379 } 380 381 #undef __FUNCT__ 382 #define __FUNCT__ "MatGetRow_SeqBAIJ" 383 int MatGetRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 384 { 385 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 386 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*idx_i,bs2,ierr; 387 MatScalar *aa,*aa_i; 388 PetscScalar *v_i; 389 390 PetscFunctionBegin; 391 bs = a->bs; 392 ai = a->i; 393 aj = a->j; 394 aa = a->a; 395 bs2 = a->bs2; 396 397 if (row < 0 || row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row out of range"); 398 399 bn = row/bs; /* Block number */ 400 bp = row % bs; /* Block Position */ 401 M = ai[bn+1] - ai[bn]; 402 *nz = bs*M; 403 404 if (v) { 405 *v = 0; 406 if (*nz) { 407 ierr = PetscMalloc((*nz)*sizeof(PetscScalar),v);CHKERRQ(ierr); 408 for (i=0; i<M; i++) { /* for each block in the block row */ 409 v_i = *v + i*bs; 410 aa_i = aa + bs2*(ai[bn] + i); 411 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 412 } 413 } 414 } 415 416 if (idx) { 417 *idx = 0; 418 if (*nz) { 419 ierr = PetscMalloc((*nz)*sizeof(int),idx);CHKERRQ(ierr); 420 for (i=0; i<M; i++) { /* for each block in the block row */ 421 idx_i = *idx + i*bs; 422 itmp = bs*aj[ai[bn] + i]; 423 for (j=0; j<bs; j++) {idx_i[j] = itmp++;} 424 } 425 } 426 } 427 PetscFunctionReturn(0); 428 } 429 430 #undef __FUNCT__ 431 #define __FUNCT__ "MatRestoreRow_SeqBAIJ" 432 int MatRestoreRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 433 { 434 int ierr; 435 436 PetscFunctionBegin; 437 if (idx) {if (*idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 438 if (v) {if (*v) {ierr = PetscFree(*v);CHKERRQ(ierr);}} 439 PetscFunctionReturn(0); 440 } 441 442 #undef __FUNCT__ 443 #define __FUNCT__ "MatTranspose_SeqBAIJ" 444 int MatTranspose_SeqBAIJ(Mat A,Mat *B) 445 { 446 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 447 Mat C; 448 int i,j,k,ierr,*aj=a->j,*ai=a->i,bs=a->bs,mbs=a->mbs,nbs=a->nbs,len,*col; 449 int *rows,*cols,bs2=a->bs2; 450 PetscScalar *array; 451 452 PetscFunctionBegin; 453 if (!B && mbs!=nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Square matrix only for in-place"); 454 ierr = PetscMalloc((1+nbs)*sizeof(int),&col);CHKERRQ(ierr); 455 ierr = PetscMemzero(col,(1+nbs)*sizeof(int));CHKERRQ(ierr); 456 457 #if defined(PETSC_USE_MAT_SINGLE) 458 ierr = PetscMalloc(a->bs2*a->nz*sizeof(PetscScalar),&array);CHKERRQ(ierr); 459 for (i=0; i<a->bs2*a->nz; i++) array[i] = (PetscScalar)a->a[i]; 460 #else 461 array = a->a; 462 #endif 463 464 for (i=0; i<ai[mbs]; i++) col[aj[i]] += 1; 465 ierr = MatCreateSeqBAIJ(A->comm,bs,A->n,A->m,PETSC_NULL,col,&C);CHKERRQ(ierr); 466 ierr = PetscFree(col);CHKERRQ(ierr); 467 ierr = PetscMalloc(2*bs*sizeof(int),&rows);CHKERRQ(ierr); 468 cols = rows + bs; 469 for (i=0; i<mbs; i++) { 470 cols[0] = i*bs; 471 for (k=1; k<bs; k++) cols[k] = cols[k-1] + 1; 472 len = ai[i+1] - ai[i]; 473 for (j=0; j<len; j++) { 474 rows[0] = (*aj++)*bs; 475 for (k=1; k<bs; k++) rows[k] = rows[k-1] + 1; 476 ierr = MatSetValues(C,bs,rows,bs,cols,array,INSERT_VALUES);CHKERRQ(ierr); 477 array += bs2; 478 } 479 } 480 ierr = PetscFree(rows);CHKERRQ(ierr); 481 #if defined(PETSC_USE_MAT_SINGLE) 482 ierr = PetscFree(array); 483 #endif 484 485 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 486 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 487 488 if (B) { 489 *B = C; 490 } else { 491 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 492 } 493 PetscFunctionReturn(0); 494 } 495 496 #undef __FUNCT__ 497 #define __FUNCT__ "MatView_SeqBAIJ_Binary" 498 static int MatView_SeqBAIJ_Binary(Mat A,PetscViewer viewer) 499 { 500 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 501 int i,fd,*col_lens,ierr,bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 502 PetscScalar *aa; 503 FILE *file; 504 505 PetscFunctionBegin; 506 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 507 ierr = PetscMalloc((4+A->m)*sizeof(int),&col_lens);CHKERRQ(ierr); 508 col_lens[0] = MAT_FILE_COOKIE; 509 510 col_lens[1] = A->m; 511 col_lens[2] = A->n; 512 col_lens[3] = a->nz*bs2; 513 514 /* store lengths of each row and write (including header) to file */ 515 count = 0; 516 for (i=0; i<a->mbs; i++) { 517 for (j=0; j<bs; j++) { 518 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 519 } 520 } 521 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,1);CHKERRQ(ierr); 522 ierr = PetscFree(col_lens);CHKERRQ(ierr); 523 524 /* store column indices (zero start index) */ 525 ierr = PetscMalloc((a->nz+1)*bs2*sizeof(int),&jj);CHKERRQ(ierr); 526 count = 0; 527 for (i=0; i<a->mbs; i++) { 528 for (j=0; j<bs; j++) { 529 for (k=a->i[i]; k<a->i[i+1]; k++) { 530 for (l=0; l<bs; l++) { 531 jj[count++] = bs*a->j[k] + l; 532 } 533 } 534 } 535 } 536 ierr = PetscBinaryWrite(fd,jj,bs2*a->nz,PETSC_INT,0);CHKERRQ(ierr); 537 ierr = PetscFree(jj);CHKERRQ(ierr); 538 539 /* store nonzero values */ 540 ierr = PetscMalloc((a->nz+1)*bs2*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 541 count = 0; 542 for (i=0; i<a->mbs; i++) { 543 for (j=0; j<bs; j++) { 544 for (k=a->i[i]; k<a->i[i+1]; k++) { 545 for (l=0; l<bs; l++) { 546 aa[count++] = a->a[bs2*k + l*bs + j]; 547 } 548 } 549 } 550 } 551 ierr = PetscBinaryWrite(fd,aa,bs2*a->nz,PETSC_SCALAR,0);CHKERRQ(ierr); 552 ierr = PetscFree(aa);CHKERRQ(ierr); 553 554 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 555 if (file) { 556 fprintf(file,"-matload_block_size %d\n",a->bs); 557 } 558 PetscFunctionReturn(0); 559 } 560 561 extern int MatMPIBAIJFactorInfo_DSCPACK(Mat,PetscViewer); 562 563 #undef __FUNCT__ 564 #define __FUNCT__ "MatView_SeqBAIJ_ASCII" 565 static int MatView_SeqBAIJ_ASCII(Mat A,PetscViewer viewer) 566 { 567 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 568 int ierr,i,j,bs = a->bs,k,l,bs2=a->bs2; 569 PetscViewerFormat format; 570 571 PetscFunctionBegin; 572 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 573 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_LONG) { 574 ierr = PetscViewerASCIIPrintf(viewer," block size is %d\n",bs);CHKERRQ(ierr); 575 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 576 Mat aij; 577 ierr = MatConvert(A,MATSEQAIJ,&aij);CHKERRQ(ierr); 578 ierr = MatView(aij,viewer);CHKERRQ(ierr); 579 ierr = MatDestroy(aij);CHKERRQ(ierr); 580 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 581 #if defined(PETSC_HAVE_DSCPACK) && !defined(PETSC_USE_SINGLE) && !defined(PETSC_USE_COMPLEX) 582 ierr = MatMPIBAIJFactorInfo_DSCPACK(A,viewer);CHKERRQ(ierr); 583 #endif 584 PetscFunctionReturn(0); 585 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 586 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 587 for (i=0; i<a->mbs; i++) { 588 for (j=0; j<bs; j++) { 589 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 590 for (k=a->i[i]; k<a->i[i+1]; k++) { 591 for (l=0; l<bs; l++) { 592 #if defined(PETSC_USE_COMPLEX) 593 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 594 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g + %gi) ",bs*a->j[k]+l, 595 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 596 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 597 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g - %gi) ",bs*a->j[k]+l, 598 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 599 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 600 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 601 } 602 #else 603 if (a->a[bs2*k + l*bs + j] != 0.0) { 604 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 605 } 606 #endif 607 } 608 } 609 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 610 } 611 } 612 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 613 } else { 614 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 615 for (i=0; i<a->mbs; i++) { 616 for (j=0; j<bs; j++) { 617 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 618 for (k=a->i[i]; k<a->i[i+1]; k++) { 619 for (l=0; l<bs; l++) { 620 #if defined(PETSC_USE_COMPLEX) 621 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 622 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g + %g i) ",bs*a->j[k]+l, 623 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 624 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 625 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g - %g i) ",bs*a->j[k]+l, 626 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 627 } else { 628 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 629 } 630 #else 631 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 632 #endif 633 } 634 } 635 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 636 } 637 } 638 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 639 } 640 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 641 PetscFunctionReturn(0); 642 } 643 644 #undef __FUNCT__ 645 #define __FUNCT__ "MatView_SeqBAIJ_Draw_Zoom" 646 static int MatView_SeqBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 647 { 648 Mat A = (Mat) Aa; 649 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 650 int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2; 651 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 652 MatScalar *aa; 653 PetscViewer viewer; 654 655 PetscFunctionBegin; 656 657 /* still need to add support for contour plot of nonzeros; see MatView_SeqAIJ_Draw_Zoom()*/ 658 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 659 660 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 661 662 /* loop over matrix elements drawing boxes */ 663 color = PETSC_DRAW_BLUE; 664 for (i=0,row=0; i<mbs; i++,row+=bs) { 665 for (j=a->i[i]; j<a->i[i+1]; j++) { 666 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 667 x_l = a->j[j]*bs; x_r = x_l + 1.0; 668 aa = a->a + j*bs2; 669 for (k=0; k<bs; k++) { 670 for (l=0; l<bs; l++) { 671 if (PetscRealPart(*aa++) >= 0.) continue; 672 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 673 } 674 } 675 } 676 } 677 color = PETSC_DRAW_CYAN; 678 for (i=0,row=0; i<mbs; i++,row+=bs) { 679 for (j=a->i[i]; j<a->i[i+1]; j++) { 680 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 681 x_l = a->j[j]*bs; x_r = x_l + 1.0; 682 aa = a->a + j*bs2; 683 for (k=0; k<bs; k++) { 684 for (l=0; l<bs; l++) { 685 if (PetscRealPart(*aa++) != 0.) continue; 686 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 687 } 688 } 689 } 690 } 691 692 color = PETSC_DRAW_RED; 693 for (i=0,row=0; i<mbs; i++,row+=bs) { 694 for (j=a->i[i]; j<a->i[i+1]; j++) { 695 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 696 x_l = a->j[j]*bs; x_r = x_l + 1.0; 697 aa = a->a + j*bs2; 698 for (k=0; k<bs; k++) { 699 for (l=0; l<bs; l++) { 700 if (PetscRealPart(*aa++) <= 0.) continue; 701 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 702 } 703 } 704 } 705 } 706 PetscFunctionReturn(0); 707 } 708 709 #undef __FUNCT__ 710 #define __FUNCT__ "MatView_SeqBAIJ_Draw" 711 static int MatView_SeqBAIJ_Draw(Mat A,PetscViewer viewer) 712 { 713 int ierr; 714 PetscReal xl,yl,xr,yr,w,h; 715 PetscDraw draw; 716 PetscTruth isnull; 717 718 PetscFunctionBegin; 719 720 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 721 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 722 723 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 724 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 725 xr += w; yr += h; xl = -w; yl = -h; 726 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 727 ierr = PetscDrawZoom(draw,MatView_SeqBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 728 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 729 PetscFunctionReturn(0); 730 } 731 732 #undef __FUNCT__ 733 #define __FUNCT__ "MatView_SeqBAIJ" 734 int MatView_SeqBAIJ(Mat A,PetscViewer viewer) 735 { 736 int ierr; 737 PetscTruth issocket,isascii,isbinary,isdraw; 738 739 PetscFunctionBegin; 740 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 741 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 742 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 743 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 744 if (issocket) { 745 SETERRQ(PETSC_ERR_SUP,"Socket viewer not supported"); 746 } else if (isascii){ 747 ierr = MatView_SeqBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 748 } else if (isbinary) { 749 ierr = MatView_SeqBAIJ_Binary(A,viewer);CHKERRQ(ierr); 750 } else if (isdraw) { 751 ierr = MatView_SeqBAIJ_Draw(A,viewer);CHKERRQ(ierr); 752 } else { 753 SETERRQ1(1,"Viewer type %s not supported by SeqBAIJ matrices",((PetscObject)viewer)->type_name); 754 } 755 PetscFunctionReturn(0); 756 } 757 758 759 #undef __FUNCT__ 760 #define __FUNCT__ "MatGetValues_SeqBAIJ" 761 int MatGetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v) 762 { 763 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 764 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 765 int *ai = a->i,*ailen = a->ilen; 766 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 767 MatScalar *ap,*aa = a->a,zero = 0.0; 768 769 PetscFunctionBegin; 770 for (k=0; k<m; k++) { /* loop over rows */ 771 row = im[k]; brow = row/bs; 772 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative row"); 773 if (row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 774 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 775 nrow = ailen[brow]; 776 for (l=0; l<n; l++) { /* loop over columns */ 777 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative column"); 778 if (in[l] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 779 col = in[l] ; 780 bcol = col/bs; 781 cidx = col%bs; 782 ridx = row%bs; 783 high = nrow; 784 low = 0; /* assume unsorted */ 785 while (high-low > 5) { 786 t = (low+high)/2; 787 if (rp[t] > bcol) high = t; 788 else low = t; 789 } 790 for (i=low; i<high; i++) { 791 if (rp[i] > bcol) break; 792 if (rp[i] == bcol) { 793 *v++ = ap[bs2*i+bs*cidx+ridx]; 794 goto finished; 795 } 796 } 797 *v++ = zero; 798 finished:; 799 } 800 } 801 PetscFunctionReturn(0); 802 } 803 804 #if defined(PETSC_USE_MAT_SINGLE) 805 #undef __FUNCT__ 806 #define __FUNCT__ "MatSetValuesBlocked_SeqBAIJ" 807 int MatSetValuesBlocked_SeqBAIJ(Mat mat,int m,int *im,int n,int *in,PetscScalar *v,InsertMode addv) 808 { 809 Mat_SeqBAIJ *b = (Mat_SeqBAIJ*)mat->data; 810 int ierr,i,N = m*n*b->bs2; 811 MatScalar *vsingle; 812 813 PetscFunctionBegin; 814 if (N > b->setvalueslen) { 815 if (b->setvaluescopy) {ierr = PetscFree(b->setvaluescopy);CHKERRQ(ierr);} 816 ierr = PetscMalloc(N*sizeof(MatScalar),&b->setvaluescopy);CHKERRQ(ierr); 817 b->setvalueslen = N; 818 } 819 vsingle = b->setvaluescopy; 820 for (i=0; i<N; i++) { 821 vsingle[i] = v[i]; 822 } 823 ierr = MatSetValuesBlocked_SeqBAIJ_MatScalar(mat,m,im,n,in,vsingle,addv);CHKERRQ(ierr); 824 PetscFunctionReturn(0); 825 } 826 #endif 827 828 829 #undef __FUNCT__ 830 #define __FUNCT__ "MatSetValuesBlocked_SeqBAIJ" 831 int MatSetValuesBlocked_SeqBAIJ_MatScalar(Mat A,int m,int *im,int n,int *in,MatScalar *v,InsertMode is) 832 { 833 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 834 int *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 835 int *imax=a->imax,*ai=a->i,*ailen=a->ilen; 836 int *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=a->bs,stepval,ierr; 837 PetscTruth roworiented=a->roworiented; 838 MatScalar *value = v,*ap,*aa = a->a,*bap; 839 840 PetscFunctionBegin; 841 if (roworiented) { 842 stepval = (n-1)*bs; 843 } else { 844 stepval = (m-1)*bs; 845 } 846 for (k=0; k<m; k++) { /* loop over added rows */ 847 row = im[k]; 848 if (row < 0) continue; 849 #if defined(PETSC_USE_BOPT_g) 850 if (row >= a->mbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 851 #endif 852 rp = aj + ai[row]; 853 ap = aa + bs2*ai[row]; 854 rmax = imax[row]; 855 nrow = ailen[row]; 856 low = 0; 857 for (l=0; l<n; l++) { /* loop over added columns */ 858 if (in[l] < 0) continue; 859 #if defined(PETSC_USE_BOPT_g) 860 if (in[l] >= a->nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 861 #endif 862 col = in[l]; 863 if (roworiented) { 864 value = v + k*(stepval+bs)*bs + l*bs; 865 } else { 866 value = v + l*(stepval+bs)*bs + k*bs; 867 } 868 if (!sorted) low = 0; high = nrow; 869 while (high-low > 7) { 870 t = (low+high)/2; 871 if (rp[t] > col) high = t; 872 else low = t; 873 } 874 for (i=low; i<high; i++) { 875 if (rp[i] > col) break; 876 if (rp[i] == col) { 877 bap = ap + bs2*i; 878 if (roworiented) { 879 if (is == ADD_VALUES) { 880 for (ii=0; ii<bs; ii++,value+=stepval) { 881 for (jj=ii; jj<bs2; jj+=bs) { 882 bap[jj] += *value++; 883 } 884 } 885 } else { 886 for (ii=0; ii<bs; ii++,value+=stepval) { 887 for (jj=ii; jj<bs2; jj+=bs) { 888 bap[jj] = *value++; 889 } 890 } 891 } 892 } else { 893 if (is == ADD_VALUES) { 894 for (ii=0; ii<bs; ii++,value+=stepval) { 895 for (jj=0; jj<bs; jj++) { 896 *bap++ += *value++; 897 } 898 } 899 } else { 900 for (ii=0; ii<bs; ii++,value+=stepval) { 901 for (jj=0; jj<bs; jj++) { 902 *bap++ = *value++; 903 } 904 } 905 } 906 } 907 goto noinsert2; 908 } 909 } 910 if (nonew == 1) goto noinsert2; 911 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 912 if (nrow >= rmax) { 913 /* there is no extra room in row, therefore enlarge */ 914 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 915 MatScalar *new_a; 916 917 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 918 919 /* malloc new storage space */ 920 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 921 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 922 new_j = (int*)(new_a + bs2*new_nz); 923 new_i = new_j + new_nz; 924 925 /* copy over old data into new slots */ 926 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 927 for (ii=row+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 928 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int));CHKERRQ(ierr); 929 len = (new_nz - CHUNKSIZE - ai[row] - nrow); 930 ierr = PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int));CHKERRQ(ierr); 931 ierr = PetscMemcpy(new_a,aa,(ai[row]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 932 ierr = PetscMemzero(new_a+bs2*(ai[row]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 933 ierr = PetscMemcpy(new_a+bs2*(ai[row]+nrow+CHUNKSIZE),aa+bs2*(ai[row]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 934 /* free up old matrix storage */ 935 ierr = PetscFree(a->a);CHKERRQ(ierr); 936 if (!a->singlemalloc) { 937 ierr = PetscFree(a->i);CHKERRQ(ierr); 938 ierr = PetscFree(a->j);CHKERRQ(ierr); 939 } 940 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 941 a->singlemalloc = PETSC_TRUE; 942 943 rp = aj + ai[row]; ap = aa + bs2*ai[row]; 944 rmax = imax[row] = imax[row] + CHUNKSIZE; 945 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 946 a->maxnz += bs2*CHUNKSIZE; 947 a->reallocs++; 948 a->nz++; 949 } 950 N = nrow++ - 1; 951 /* shift up all the later entries in this row */ 952 for (ii=N; ii>=i; ii--) { 953 rp[ii+1] = rp[ii]; 954 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 955 } 956 if (N >= i) { 957 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 958 } 959 rp[i] = col; 960 bap = ap + bs2*i; 961 if (roworiented) { 962 for (ii=0; ii<bs; ii++,value+=stepval) { 963 for (jj=ii; jj<bs2; jj+=bs) { 964 bap[jj] = *value++; 965 } 966 } 967 } else { 968 for (ii=0; ii<bs; ii++,value+=stepval) { 969 for (jj=0; jj<bs; jj++) { 970 *bap++ = *value++; 971 } 972 } 973 } 974 noinsert2:; 975 low = i; 976 } 977 ailen[row] = nrow; 978 } 979 PetscFunctionReturn(0); 980 } 981 982 #undef __FUNCT__ 983 #define __FUNCT__ "MatAssemblyEnd_SeqBAIJ" 984 int MatAssemblyEnd_SeqBAIJ(Mat A,MatAssemblyType mode) 985 { 986 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 987 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 988 int m = A->m,*ip,N,*ailen = a->ilen; 989 int mbs = a->mbs,bs2 = a->bs2,rmax = 0,ierr; 990 MatScalar *aa = a->a,*ap; 991 #if defined(PETSC_HAVE_DSCPACK) 992 PetscTruth flag; 993 #endif 994 995 PetscFunctionBegin; 996 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 997 998 if (m) rmax = ailen[0]; 999 for (i=1; i<mbs; i++) { 1000 /* move each row back by the amount of empty slots (fshift) before it*/ 1001 fshift += imax[i-1] - ailen[i-1]; 1002 rmax = PetscMax(rmax,ailen[i]); 1003 if (fshift) { 1004 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 1005 N = ailen[i]; 1006 for (j=0; j<N; j++) { 1007 ip[j-fshift] = ip[j]; 1008 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1009 } 1010 } 1011 ai[i] = ai[i-1] + ailen[i-1]; 1012 } 1013 if (mbs) { 1014 fshift += imax[mbs-1] - ailen[mbs-1]; 1015 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 1016 } 1017 /* reset ilen and imax for each row */ 1018 for (i=0; i<mbs; i++) { 1019 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1020 } 1021 a->nz = ai[mbs]; 1022 1023 /* diagonals may have moved, so kill the diagonal pointers */ 1024 if (fshift && a->diag) { 1025 ierr = PetscFree(a->diag);CHKERRQ(ierr); 1026 PetscLogObjectMemory(A,-(mbs+1)*sizeof(int)); 1027 a->diag = 0; 1028 } 1029 PetscLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n",m,A->n,a->bs,fshift*bs2,a->nz*bs2); 1030 PetscLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Number of mallocs during MatSetValues is %d\n",a->reallocs); 1031 PetscLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 1032 a->reallocs = 0; 1033 A->info.nz_unneeded = (PetscReal)fshift*bs2; 1034 1035 #if defined(PETSC_HAVE_DSCPACK) 1036 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_baij_dscpack",&flag);CHKERRQ(ierr); 1037 if (flag) { ierr = MatUseDSCPACK_MPIBAIJ(A);CHKERRQ(ierr); } 1038 #endif 1039 1040 PetscFunctionReturn(0); 1041 } 1042 1043 1044 1045 /* 1046 This function returns an array of flags which indicate the locations of contiguous 1047 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 1048 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 1049 Assume: sizes should be long enough to hold all the values. 1050 */ 1051 #undef __FUNCT__ 1052 #define __FUNCT__ "MatZeroRows_SeqBAIJ_Check_Blocks" 1053 static int MatZeroRows_SeqBAIJ_Check_Blocks(int idx[],int n,int bs,int sizes[], int *bs_max) 1054 { 1055 int i,j,k,row; 1056 PetscTruth flg; 1057 1058 PetscFunctionBegin; 1059 for (i=0,j=0; i<n; j++) { 1060 row = idx[i]; 1061 if (row%bs!=0) { /* Not the begining of a block */ 1062 sizes[j] = 1; 1063 i++; 1064 } else if (i+bs > n) { /* complete block doesn't exist (at idx end) */ 1065 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 1066 i++; 1067 } else { /* Begining of the block, so check if the complete block exists */ 1068 flg = PETSC_TRUE; 1069 for (k=1; k<bs; k++) { 1070 if (row+k != idx[i+k]) { /* break in the block */ 1071 flg = PETSC_FALSE; 1072 break; 1073 } 1074 } 1075 if (flg == PETSC_TRUE) { /* No break in the bs */ 1076 sizes[j] = bs; 1077 i+= bs; 1078 } else { 1079 sizes[j] = 1; 1080 i++; 1081 } 1082 } 1083 } 1084 *bs_max = j; 1085 PetscFunctionReturn(0); 1086 } 1087 1088 #undef __FUNCT__ 1089 #define __FUNCT__ "MatZeroRows_SeqBAIJ" 1090 int MatZeroRows_SeqBAIJ(Mat A,IS is,PetscScalar *diag) 1091 { 1092 Mat_SeqBAIJ *baij=(Mat_SeqBAIJ*)A->data; 1093 int ierr,i,j,k,count,is_n,*is_idx,*rows; 1094 int bs=baij->bs,bs2=baij->bs2,*sizes,row,bs_max; 1095 PetscScalar zero = 0.0; 1096 MatScalar *aa; 1097 1098 PetscFunctionBegin; 1099 /* Make a copy of the IS and sort it */ 1100 ierr = ISGetLocalSize(is,&is_n);CHKERRQ(ierr); 1101 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 1102 1103 /* allocate memory for rows,sizes */ 1104 ierr = PetscMalloc((3*is_n+1)*sizeof(int),&rows);CHKERRQ(ierr); 1105 sizes = rows + is_n; 1106 1107 /* copy IS values to rows, and sort them */ 1108 for (i=0; i<is_n; i++) { rows[i] = is_idx[i]; } 1109 ierr = PetscSortInt(is_n,rows);CHKERRQ(ierr); 1110 if (baij->keepzeroedrows) { 1111 for (i=0; i<is_n; i++) { sizes[i] = 1; } 1112 bs_max = is_n; 1113 } else { 1114 ierr = MatZeroRows_SeqBAIJ_Check_Blocks(rows,is_n,bs,sizes,&bs_max);CHKERRQ(ierr); 1115 } 1116 ierr = ISRestoreIndices(is,&is_idx);CHKERRQ(ierr); 1117 1118 for (i=0,j=0; i<bs_max; j+=sizes[i],i++) { 1119 row = rows[j]; 1120 if (row < 0 || row > A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %d out of range",row); 1121 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1122 aa = baij->a + baij->i[row/bs]*bs2 + (row%bs); 1123 if (sizes[i] == bs && !baij->keepzeroedrows) { 1124 if (diag) { 1125 if (baij->ilen[row/bs] > 0) { 1126 baij->ilen[row/bs] = 1; 1127 baij->j[baij->i[row/bs]] = row/bs; 1128 ierr = PetscMemzero(aa,count*bs*sizeof(MatScalar));CHKERRQ(ierr); 1129 } 1130 /* Now insert all the diagonal values for this bs */ 1131 for (k=0; k<bs; k++) { 1132 ierr = (*A->ops->setvalues)(A,1,rows+j+k,1,rows+j+k,diag,INSERT_VALUES);CHKERRQ(ierr); 1133 } 1134 } else { /* (!diag) */ 1135 baij->ilen[row/bs] = 0; 1136 } /* end (!diag) */ 1137 } else { /* (sizes[i] != bs) */ 1138 #if defined (PETSC_USE_DEBUG) 1139 if (sizes[i] != 1) SETERRQ(1,"Internal Error. Value should be 1"); 1140 #endif 1141 for (k=0; k<count; k++) { 1142 aa[0] = zero; 1143 aa += bs; 1144 } 1145 if (diag) { 1146 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 1147 } 1148 } 1149 } 1150 1151 ierr = PetscFree(rows);CHKERRQ(ierr); 1152 ierr = MatAssemblyEnd_SeqBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1153 PetscFunctionReturn(0); 1154 } 1155 1156 #undef __FUNCT__ 1157 #define __FUNCT__ "MatSetValues_SeqBAIJ" 1158 int MatSetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v,InsertMode is) 1159 { 1160 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1161 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 1162 int *imax=a->imax,*ai=a->i,*ailen=a->ilen; 1163 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 1164 int ridx,cidx,bs2=a->bs2,ierr; 1165 PetscTruth roworiented=a->roworiented; 1166 MatScalar *ap,value,*aa=a->a,*bap; 1167 1168 PetscFunctionBegin; 1169 for (k=0; k<m; k++) { /* loop over added rows */ 1170 row = im[k]; brow = row/bs; 1171 if (row < 0) continue; 1172 #if defined(PETSC_USE_BOPT_g) 1173 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m); 1174 #endif 1175 rp = aj + ai[brow]; 1176 ap = aa + bs2*ai[brow]; 1177 rmax = imax[brow]; 1178 nrow = ailen[brow]; 1179 low = 0; 1180 for (l=0; l<n; l++) { /* loop over added columns */ 1181 if (in[l] < 0) continue; 1182 #if defined(PETSC_USE_BOPT_g) 1183 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->n); 1184 #endif 1185 col = in[l]; bcol = col/bs; 1186 ridx = row % bs; cidx = col % bs; 1187 if (roworiented) { 1188 value = v[l + k*n]; 1189 } else { 1190 value = v[k + l*m]; 1191 } 1192 if (!sorted) low = 0; high = nrow; 1193 while (high-low > 7) { 1194 t = (low+high)/2; 1195 if (rp[t] > bcol) high = t; 1196 else low = t; 1197 } 1198 for (i=low; i<high; i++) { 1199 if (rp[i] > bcol) break; 1200 if (rp[i] == bcol) { 1201 bap = ap + bs2*i + bs*cidx + ridx; 1202 if (is == ADD_VALUES) *bap += value; 1203 else *bap = value; 1204 goto noinsert1; 1205 } 1206 } 1207 if (nonew == 1) goto noinsert1; 1208 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 1209 if (nrow >= rmax) { 1210 /* there is no extra room in row, therefore enlarge */ 1211 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 1212 MatScalar *new_a; 1213 1214 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 1215 1216 /* Malloc new storage space */ 1217 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 1218 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 1219 new_j = (int*)(new_a + bs2*new_nz); 1220 new_i = new_j + new_nz; 1221 1222 /* copy over old data into new slots */ 1223 for (ii=0; ii<brow+1; ii++) {new_i[ii] = ai[ii];} 1224 for (ii=brow+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 1225 ierr = PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int));CHKERRQ(ierr); 1226 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 1227 ierr = PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow,len*sizeof(int));CHKERRQ(ierr); 1228 ierr = PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1229 ierr = PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 1230 ierr = PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE),aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 1231 /* free up old matrix storage */ 1232 ierr = PetscFree(a->a);CHKERRQ(ierr); 1233 if (!a->singlemalloc) { 1234 ierr = PetscFree(a->i);CHKERRQ(ierr); 1235 ierr = PetscFree(a->j);CHKERRQ(ierr); 1236 } 1237 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 1238 a->singlemalloc = PETSC_TRUE; 1239 1240 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 1241 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 1242 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 1243 a->maxnz += bs2*CHUNKSIZE; 1244 a->reallocs++; 1245 a->nz++; 1246 } 1247 N = nrow++ - 1; 1248 /* shift up all the later entries in this row */ 1249 for (ii=N; ii>=i; ii--) { 1250 rp[ii+1] = rp[ii]; 1251 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 1252 } 1253 if (N>=i) { 1254 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1255 } 1256 rp[i] = bcol; 1257 ap[bs2*i + bs*cidx + ridx] = value; 1258 noinsert1:; 1259 low = i; 1260 } 1261 ailen[brow] = nrow; 1262 } 1263 PetscFunctionReturn(0); 1264 } 1265 1266 1267 #undef __FUNCT__ 1268 #define __FUNCT__ "MatILUFactor_SeqBAIJ" 1269 int MatILUFactor_SeqBAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1270 { 1271 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)inA->data; 1272 Mat outA; 1273 int ierr; 1274 PetscTruth row_identity,col_identity; 1275 1276 PetscFunctionBegin; 1277 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels = 0 supported for in-place ILU"); 1278 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1279 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1280 if (!row_identity || !col_identity) { 1281 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1282 } 1283 1284 outA = inA; 1285 inA->factor = FACTOR_LU; 1286 1287 if (!a->diag) { 1288 ierr = MatMarkDiagonal_SeqBAIJ(inA);CHKERRQ(ierr); 1289 } 1290 1291 a->row = row; 1292 a->col = col; 1293 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1294 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1295 1296 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1297 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1298 PetscLogObjectParent(inA,a->icol); 1299 1300 /* 1301 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 1302 for ILU(0) factorization with natural ordering 1303 */ 1304 if (a->bs < 8) { 1305 ierr = MatSeqBAIJ_UpdateFactorNumeric_NaturalOrdering(inA);CHKERRQ(ierr); 1306 } else { 1307 if (!a->solve_work) { 1308 ierr = PetscMalloc((inA->m+a->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1309 PetscLogObjectMemory(inA,(inA->m+a->bs)*sizeof(PetscScalar)); 1310 } 1311 } 1312 1313 ierr = MatLUFactorNumeric(inA,&outA);CHKERRQ(ierr); 1314 1315 PetscFunctionReturn(0); 1316 } 1317 #undef __FUNCT__ 1318 #define __FUNCT__ "MatPrintHelp_SeqBAIJ" 1319 int MatPrintHelp_SeqBAIJ(Mat A) 1320 { 1321 static PetscTruth called = PETSC_FALSE; 1322 MPI_Comm comm = A->comm; 1323 int ierr; 1324 1325 PetscFunctionBegin; 1326 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1327 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQBAIJ and MATMPIBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1328 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1329 PetscFunctionReturn(0); 1330 } 1331 1332 EXTERN_C_BEGIN 1333 #undef __FUNCT__ 1334 #define __FUNCT__ "MatSeqBAIJSetColumnIndices_SeqBAIJ" 1335 int MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat,int *indices) 1336 { 1337 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1338 int i,nz,nbs; 1339 1340 PetscFunctionBegin; 1341 nz = baij->maxnz/baij->bs2; 1342 nbs = baij->nbs; 1343 for (i=0; i<nz; i++) { 1344 baij->j[i] = indices[i]; 1345 } 1346 baij->nz = nz; 1347 for (i=0; i<nbs; i++) { 1348 baij->ilen[i] = baij->imax[i]; 1349 } 1350 1351 PetscFunctionReturn(0); 1352 } 1353 EXTERN_C_END 1354 1355 #undef __FUNCT__ 1356 #define __FUNCT__ "MatSeqBAIJSetColumnIndices" 1357 /*@ 1358 MatSeqBAIJSetColumnIndices - Set the column indices for all the rows 1359 in the matrix. 1360 1361 Input Parameters: 1362 + mat - the SeqBAIJ matrix 1363 - indices - the column indices 1364 1365 Level: advanced 1366 1367 Notes: 1368 This can be called if you have precomputed the nonzero structure of the 1369 matrix and want to provide it to the matrix object to improve the performance 1370 of the MatSetValues() operation. 1371 1372 You MUST have set the correct numbers of nonzeros per row in the call to 1373 MatCreateSeqBAIJ(). 1374 1375 MUST be called before any calls to MatSetValues(); 1376 1377 @*/ 1378 int MatSeqBAIJSetColumnIndices(Mat mat,int *indices) 1379 { 1380 int ierr,(*f)(Mat,int *); 1381 1382 PetscFunctionBegin; 1383 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1384 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1385 if (f) { 1386 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1387 } else { 1388 SETERRQ(1,"Wrong type of matrix to set column indices"); 1389 } 1390 PetscFunctionReturn(0); 1391 } 1392 1393 #undef __FUNCT__ 1394 #define __FUNCT__ "MatGetRowMax_SeqBAIJ" 1395 int MatGetRowMax_SeqBAIJ(Mat A,Vec v) 1396 { 1397 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1398 int ierr,i,j,n,row,bs,*ai,*aj,mbs; 1399 PetscReal atmp; 1400 PetscScalar *x,zero = 0.0; 1401 MatScalar *aa; 1402 int ncols,brow,krow,kcol; 1403 1404 PetscFunctionBegin; 1405 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1406 bs = a->bs; 1407 aa = a->a; 1408 ai = a->i; 1409 aj = a->j; 1410 mbs = a->mbs; 1411 1412 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1413 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1414 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1415 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1416 for (i=0; i<mbs; i++) { 1417 ncols = ai[1] - ai[0]; ai++; 1418 brow = bs*i; 1419 for (j=0; j<ncols; j++){ 1420 /* bcol = bs*(*aj); */ 1421 for (kcol=0; kcol<bs; kcol++){ 1422 for (krow=0; krow<bs; krow++){ 1423 atmp = PetscAbsScalar(*aa); aa++; 1424 row = brow + krow; /* row index */ 1425 /* printf("val[%d,%d]: %g\n",row,bcol+kcol,atmp); */ 1426 if (PetscAbsScalar(x[row]) < atmp) x[row] = atmp; 1427 } 1428 } 1429 aj++; 1430 } 1431 } 1432 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1433 PetscFunctionReturn(0); 1434 } 1435 1436 #undef __FUNCT__ 1437 #define __FUNCT__ "MatSetUpPreallocation_SeqBAIJ" 1438 int MatSetUpPreallocation_SeqBAIJ(Mat A) 1439 { 1440 int ierr; 1441 1442 PetscFunctionBegin; 1443 ierr = MatSeqBAIJSetPreallocation(A,1,PETSC_DEFAULT,0);CHKERRQ(ierr); 1444 PetscFunctionReturn(0); 1445 } 1446 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "MatGetArray_SeqBAIJ" 1449 int MatGetArray_SeqBAIJ(Mat A,PetscScalar **array) 1450 { 1451 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1452 PetscFunctionBegin; 1453 *array = a->a; 1454 PetscFunctionReturn(0); 1455 } 1456 1457 #undef __FUNCT__ 1458 #define __FUNCT__ "MatRestoreArray_SeqBAIJ" 1459 int MatRestoreArray_SeqBAIJ(Mat A,PetscScalar **array) 1460 { 1461 PetscFunctionBegin; 1462 PetscFunctionReturn(0); 1463 } 1464 1465 /* -------------------------------------------------------------------*/ 1466 static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ, 1467 MatGetRow_SeqBAIJ, 1468 MatRestoreRow_SeqBAIJ, 1469 MatMult_SeqBAIJ_N, 1470 MatMultAdd_SeqBAIJ_N, 1471 MatMultTranspose_SeqBAIJ, 1472 MatMultTransposeAdd_SeqBAIJ, 1473 MatSolve_SeqBAIJ_N, 1474 0, 1475 0, 1476 0, 1477 MatLUFactor_SeqBAIJ, 1478 0, 1479 0, 1480 MatTranspose_SeqBAIJ, 1481 MatGetInfo_SeqBAIJ, 1482 MatEqual_SeqBAIJ, 1483 MatGetDiagonal_SeqBAIJ, 1484 MatDiagonalScale_SeqBAIJ, 1485 MatNorm_SeqBAIJ, 1486 0, 1487 MatAssemblyEnd_SeqBAIJ, 1488 0, 1489 MatSetOption_SeqBAIJ, 1490 MatZeroEntries_SeqBAIJ, 1491 MatZeroRows_SeqBAIJ, 1492 MatLUFactorSymbolic_SeqBAIJ, 1493 MatLUFactorNumeric_SeqBAIJ_N, 1494 0, 1495 0, 1496 MatSetUpPreallocation_SeqBAIJ, 1497 MatILUFactorSymbolic_SeqBAIJ, 1498 0, 1499 MatGetArray_SeqBAIJ, 1500 MatRestoreArray_SeqBAIJ, 1501 MatDuplicate_SeqBAIJ, 1502 0, 1503 0, 1504 MatILUFactor_SeqBAIJ, 1505 0, 1506 0, 1507 MatGetSubMatrices_SeqBAIJ, 1508 MatIncreaseOverlap_SeqBAIJ, 1509 MatGetValues_SeqBAIJ, 1510 0, 1511 MatPrintHelp_SeqBAIJ, 1512 MatScale_SeqBAIJ, 1513 0, 1514 0, 1515 0, 1516 MatGetBlockSize_SeqBAIJ, 1517 MatGetRowIJ_SeqBAIJ, 1518 MatRestoreRowIJ_SeqBAIJ, 1519 0, 1520 0, 1521 0, 1522 0, 1523 0, 1524 0, 1525 MatSetValuesBlocked_SeqBAIJ, 1526 MatGetSubMatrix_SeqBAIJ, 1527 MatDestroy_SeqBAIJ, 1528 MatView_SeqBAIJ, 1529 MatGetPetscMaps_Petsc, 1530 0, 1531 0, 1532 0, 1533 0, 1534 0, 1535 0, 1536 MatGetRowMax_SeqBAIJ, 1537 MatConvert_Basic}; 1538 1539 EXTERN_C_BEGIN 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "MatStoreValues_SeqBAIJ" 1542 int MatStoreValues_SeqBAIJ(Mat mat) 1543 { 1544 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1545 int nz = aij->i[mat->m]*aij->bs*aij->bs2; 1546 int ierr; 1547 1548 PetscFunctionBegin; 1549 if (aij->nonew != 1) { 1550 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1551 } 1552 1553 /* allocate space for values if not already there */ 1554 if (!aij->saved_values) { 1555 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1556 } 1557 1558 /* copy values over */ 1559 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1560 PetscFunctionReturn(0); 1561 } 1562 EXTERN_C_END 1563 1564 EXTERN_C_BEGIN 1565 #undef __FUNCT__ 1566 #define __FUNCT__ "MatRetrieveValues_SeqBAIJ" 1567 int MatRetrieveValues_SeqBAIJ(Mat mat) 1568 { 1569 Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data; 1570 int nz = aij->i[mat->m]*aij->bs*aij->bs2,ierr; 1571 1572 PetscFunctionBegin; 1573 if (aij->nonew != 1) { 1574 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1575 } 1576 if (!aij->saved_values) { 1577 SETERRQ(1,"Must call MatStoreValues(A);first"); 1578 } 1579 1580 /* copy values over */ 1581 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1582 PetscFunctionReturn(0); 1583 } 1584 EXTERN_C_END 1585 1586 EXTERN_C_BEGIN 1587 extern int MatConvert_SeqBAIJ_SeqAIJ(Mat,MatType,Mat*); 1588 EXTERN_C_END 1589 1590 EXTERN_C_BEGIN 1591 #undef __FUNCT__ 1592 #define __FUNCT__ "MatCreate_SeqBAIJ" 1593 int MatCreate_SeqBAIJ(Mat B) 1594 { 1595 int ierr,size; 1596 Mat_SeqBAIJ *b; 1597 1598 PetscFunctionBegin; 1599 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1600 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1601 1602 B->m = B->M = PetscMax(B->m,B->M); 1603 B->n = B->N = PetscMax(B->n,B->N); 1604 ierr = PetscNew(Mat_SeqBAIJ,&b);CHKERRQ(ierr); 1605 B->data = (void*)b; 1606 ierr = PetscMemzero(b,sizeof(Mat_SeqBAIJ));CHKERRQ(ierr); 1607 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1608 B->factor = 0; 1609 B->lupivotthreshold = 1.0; 1610 B->mapping = 0; 1611 b->row = 0; 1612 b->col = 0; 1613 b->icol = 0; 1614 b->reallocs = 0; 1615 b->saved_values = 0; 1616 #if defined(PETSC_USE_MAT_SINGLE) 1617 b->setvalueslen = 0; 1618 b->setvaluescopy = PETSC_NULL; 1619 #endif 1620 b->single_precision_solves = PETSC_FALSE; 1621 1622 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1623 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1624 1625 b->sorted = PETSC_FALSE; 1626 b->roworiented = PETSC_TRUE; 1627 b->nonew = 0; 1628 b->diag = 0; 1629 b->solve_work = 0; 1630 b->mult_work = 0; 1631 B->spptr = 0; 1632 B->info.nz_unneeded = (PetscReal)b->maxnz; 1633 b->keepzeroedrows = PETSC_FALSE; 1634 1635 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1636 "MatStoreValues_SeqBAIJ", 1637 MatStoreValues_SeqBAIJ);CHKERRQ(ierr); 1638 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1639 "MatRetrieveValues_SeqBAIJ", 1640 MatRetrieveValues_SeqBAIJ);CHKERRQ(ierr); 1641 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqBAIJSetColumnIndices_C", 1642 "MatSeqBAIJSetColumnIndices_SeqBAIJ", 1643 MatSeqBAIJSetColumnIndices_SeqBAIJ);CHKERRQ(ierr); 1644 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_SeqBAIJ_SeqAIJ_C", 1645 "MatConvert_SeqBAIJ_SeqAIJ", 1646 MatConvert_SeqBAIJ_SeqAIJ);CHKERRQ(ierr); 1647 PetscFunctionReturn(0); 1648 } 1649 EXTERN_C_END 1650 1651 #undef __FUNCT__ 1652 #define __FUNCT__ "MatDuplicate_SeqBAIJ" 1653 int MatDuplicate_SeqBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1654 { 1655 Mat C; 1656 Mat_SeqBAIJ *c,*a = (Mat_SeqBAIJ*)A->data; 1657 int i,len,mbs = a->mbs,nz = a->nz,bs2 =a->bs2,ierr; 1658 1659 PetscFunctionBegin; 1660 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1661 1662 *B = 0; 1663 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 1664 ierr = MatSetType(C,MATSEQBAIJ);CHKERRQ(ierr); 1665 c = (Mat_SeqBAIJ*)C->data; 1666 1667 c->bs = a->bs; 1668 c->bs2 = a->bs2; 1669 c->mbs = a->mbs; 1670 c->nbs = a->nbs; 1671 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1672 1673 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 1674 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 1675 for (i=0; i<mbs; i++) { 1676 c->imax[i] = a->imax[i]; 1677 c->ilen[i] = a->ilen[i]; 1678 } 1679 1680 /* allocate the matrix space */ 1681 c->singlemalloc = PETSC_TRUE; 1682 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1683 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 1684 c->j = (int*)(c->a + nz*bs2); 1685 c->i = c->j + nz; 1686 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1687 if (mbs > 0) { 1688 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1689 if (cpvalues == MAT_COPY_VALUES) { 1690 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1691 } else { 1692 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1693 } 1694 } 1695 1696 PetscLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1697 c->sorted = a->sorted; 1698 c->roworiented = a->roworiented; 1699 c->nonew = a->nonew; 1700 1701 if (a->diag) { 1702 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 1703 PetscLogObjectMemory(C,(mbs+1)*sizeof(int)); 1704 for (i=0; i<mbs; i++) { 1705 c->diag[i] = a->diag[i]; 1706 } 1707 } else c->diag = 0; 1708 c->nz = a->nz; 1709 c->maxnz = a->maxnz; 1710 c->solve_work = 0; 1711 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1712 c->mult_work = 0; 1713 C->preallocated = PETSC_TRUE; 1714 C->assembled = PETSC_TRUE; 1715 *B = C; 1716 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1717 PetscFunctionReturn(0); 1718 } 1719 1720 EXTERN_C_BEGIN 1721 #undef __FUNCT__ 1722 #define __FUNCT__ "MatLoad_SeqBAIJ" 1723 int MatLoad_SeqBAIJ(PetscViewer viewer,MatType type,Mat *A) 1724 { 1725 Mat_SeqBAIJ *a; 1726 Mat B; 1727 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1728 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1729 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1730 int *masked,nmask,tmp,bs2,ishift; 1731 PetscScalar *aa; 1732 MPI_Comm comm = ((PetscObject)viewer)->comm; 1733 1734 PetscFunctionBegin; 1735 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1736 bs2 = bs*bs; 1737 1738 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1739 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1740 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1741 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1742 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1743 M = header[1]; N = header[2]; nz = header[3]; 1744 1745 if (header[3] < 0) { 1746 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqBAIJ"); 1747 } 1748 1749 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1750 1751 /* 1752 This code adds extra rows to make sure the number of rows is 1753 divisible by the blocksize 1754 */ 1755 mbs = M/bs; 1756 extra_rows = bs - M + bs*(mbs); 1757 if (extra_rows == bs) extra_rows = 0; 1758 else mbs++; 1759 if (extra_rows) { 1760 PetscLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1761 } 1762 1763 /* read in row lengths */ 1764 ierr = PetscMalloc((M+extra_rows)*sizeof(int),&rowlengths);CHKERRQ(ierr); 1765 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1766 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1767 1768 /* read in column indices */ 1769 ierr = PetscMalloc((nz+extra_rows)*sizeof(int),&jj);CHKERRQ(ierr); 1770 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1771 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1772 1773 /* loop over row lengths determining block row lengths */ 1774 ierr = PetscMalloc(mbs*sizeof(int),&browlengths);CHKERRQ(ierr); 1775 ierr = PetscMemzero(browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1776 ierr = PetscMalloc(2*mbs*sizeof(int),&mask);CHKERRQ(ierr); 1777 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1778 masked = mask + mbs; 1779 rowcount = 0; nzcount = 0; 1780 for (i=0; i<mbs; i++) { 1781 nmask = 0; 1782 for (j=0; j<bs; j++) { 1783 kmax = rowlengths[rowcount]; 1784 for (k=0; k<kmax; k++) { 1785 tmp = jj[nzcount++]/bs; 1786 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1787 } 1788 rowcount++; 1789 } 1790 browlengths[i] += nmask; 1791 /* zero out the mask elements we set */ 1792 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1793 } 1794 1795 /* create our matrix */ 1796 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1797 B = *A; 1798 a = (Mat_SeqBAIJ*)B->data; 1799 1800 /* set matrix "i" values */ 1801 a->i[0] = 0; 1802 for (i=1; i<= mbs; i++) { 1803 a->i[i] = a->i[i-1] + browlengths[i-1]; 1804 a->ilen[i-1] = browlengths[i-1]; 1805 } 1806 a->nz = 0; 1807 for (i=0; i<mbs; i++) a->nz += browlengths[i]; 1808 1809 /* read in nonzero values */ 1810 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 1811 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1812 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1813 1814 /* set "a" and "j" values into matrix */ 1815 nzcount = 0; jcount = 0; 1816 for (i=0; i<mbs; i++) { 1817 nzcountb = nzcount; 1818 nmask = 0; 1819 for (j=0; j<bs; j++) { 1820 kmax = rowlengths[i*bs+j]; 1821 for (k=0; k<kmax; k++) { 1822 tmp = jj[nzcount++]/bs; 1823 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1824 } 1825 } 1826 /* sort the masked values */ 1827 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1828 1829 /* set "j" values into matrix */ 1830 maskcount = 1; 1831 for (j=0; j<nmask; j++) { 1832 a->j[jcount++] = masked[j]; 1833 mask[masked[j]] = maskcount++; 1834 } 1835 /* set "a" values into matrix */ 1836 ishift = bs2*a->i[i]; 1837 for (j=0; j<bs; j++) { 1838 kmax = rowlengths[i*bs+j]; 1839 for (k=0; k<kmax; k++) { 1840 tmp = jj[nzcountb]/bs ; 1841 block = mask[tmp] - 1; 1842 point = jj[nzcountb] - bs*tmp; 1843 idx = ishift + bs2*block + j + bs*point; 1844 a->a[idx] = (MatScalar)aa[nzcountb++]; 1845 } 1846 } 1847 /* zero out the mask elements we set */ 1848 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1849 } 1850 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1851 1852 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1853 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1854 ierr = PetscFree(aa);CHKERRQ(ierr); 1855 ierr = PetscFree(jj);CHKERRQ(ierr); 1856 ierr = PetscFree(mask);CHKERRQ(ierr); 1857 1858 B->assembled = PETSC_TRUE; 1859 1860 ierr = MatView_Private(B);CHKERRQ(ierr); 1861 PetscFunctionReturn(0); 1862 } 1863 EXTERN_C_END 1864 1865 #undef __FUNCT__ 1866 #define __FUNCT__ "MatCreateSeqBAIJ" 1867 /*@C 1868 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1869 compressed row) format. For good matrix assembly performance the 1870 user should preallocate the matrix storage by setting the parameter nz 1871 (or the array nnz). By setting these parameters accurately, performance 1872 during matrix assembly can be increased by more than a factor of 50. 1873 1874 Collective on MPI_Comm 1875 1876 Input Parameters: 1877 + comm - MPI communicator, set to PETSC_COMM_SELF 1878 . bs - size of block 1879 . m - number of rows 1880 . n - number of columns 1881 . nz - number of nonzero blocks per block row (same for all rows) 1882 - nnz - array containing the number of nonzero blocks in the various block rows 1883 (possibly different for each block row) or PETSC_NULL 1884 1885 Output Parameter: 1886 . A - the matrix 1887 1888 Options Database Keys: 1889 . -mat_no_unroll - uses code that does not unroll the loops in the 1890 block calculations (much slower) 1891 . -mat_block_size - size of the blocks to use 1892 1893 Level: intermediate 1894 1895 Notes: 1896 A nonzero block is any block that as 1 or more nonzeros in it 1897 1898 The block AIJ format is fully compatible with standard Fortran 77 1899 storage. That is, the stored row and column indices can begin at 1900 either one (as in Fortran) or zero. See the users' manual for details. 1901 1902 Specify the preallocated storage with either nz or nnz (not both). 1903 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1904 allocation. For additional details, see the users manual chapter on 1905 matrices. 1906 1907 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1908 @*/ 1909 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1910 { 1911 int ierr; 1912 1913 PetscFunctionBegin; 1914 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1915 ierr = MatSetType(*A,MATSEQBAIJ);CHKERRQ(ierr); 1916 ierr = MatSeqBAIJSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1917 PetscFunctionReturn(0); 1918 } 1919 1920 #undef __FUNCT__ 1921 #define __FUNCT__ "MatSeqBAIJSetPreallocation" 1922 /*@C 1923 MatSeqBAIJSetPreallocation - Sets the block size and expected nonzeros 1924 per row in the matrix. For good matrix assembly performance the 1925 user should preallocate the matrix storage by setting the parameter nz 1926 (or the array nnz). By setting these parameters accurately, performance 1927 during matrix assembly can be increased by more than a factor of 50. 1928 1929 Collective on MPI_Comm 1930 1931 Input Parameters: 1932 + A - the matrix 1933 . bs - size of block 1934 . nz - number of block nonzeros per block row (same for all rows) 1935 - nnz - array containing the number of block nonzeros in the various block rows 1936 (possibly different for each block row) or PETSC_NULL 1937 1938 Options Database Keys: 1939 . -mat_no_unroll - uses code that does not unroll the loops in the 1940 block calculations (much slower) 1941 . -mat_block_size - size of the blocks to use 1942 1943 Level: intermediate 1944 1945 Notes: 1946 The block AIJ format is fully compatible with standard Fortran 77 1947 storage. That is, the stored row and column indices can begin at 1948 either one (as in Fortran) or zero. See the users' manual for details. 1949 1950 Specify the preallocated storage with either nz or nnz (not both). 1951 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1952 allocation. For additional details, see the users manual chapter on 1953 matrices. 1954 1955 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1956 @*/ 1957 int MatSeqBAIJSetPreallocation(Mat B,int bs,int nz,int *nnz) 1958 { 1959 Mat_SeqBAIJ *b; 1960 int i,len,ierr,mbs,nbs,bs2,newbs = bs; 1961 PetscTruth flg; 1962 1963 PetscFunctionBegin; 1964 ierr = PetscTypeCompare((PetscObject)B,MATSEQBAIJ,&flg);CHKERRQ(ierr); 1965 if (!flg) PetscFunctionReturn(0); 1966 1967 B->preallocated = PETSC_TRUE; 1968 ierr = PetscOptionsGetInt(B->prefix,"-mat_block_size",&newbs,PETSC_NULL);CHKERRQ(ierr); 1969 if (nnz && newbs != bs) { 1970 SETERRQ(1,"Cannot change blocksize from command line if setting nnz"); 1971 } 1972 bs = newbs; 1973 1974 mbs = B->m/bs; 1975 nbs = B->n/bs; 1976 bs2 = bs*bs; 1977 1978 if (mbs*bs!=B->m || nbs*bs!=B->n) { 1979 SETERRQ3(PETSC_ERR_ARG_SIZ,"Number rows %d, cols %d must be divisible by blocksize %d",B->m,B->n,bs); 1980 } 1981 1982 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1983 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 1984 if (nnz) { 1985 for (i=0; i<mbs; i++) { 1986 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1987 if (nnz[i] > nbs) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %d value %d rowlength %d",i,nnz[i],nbs); 1988 } 1989 } 1990 1991 b = (Mat_SeqBAIJ*)B->data; 1992 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1993 B->ops->solve = MatSolve_SeqBAIJ_Update; 1994 B->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_Update; 1995 if (!flg) { 1996 switch (bs) { 1997 case 1: 1998 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1999 B->ops->mult = MatMult_SeqBAIJ_1; 2000 B->ops->multadd = MatMultAdd_SeqBAIJ_1; 2001 break; 2002 case 2: 2003 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 2004 B->ops->mult = MatMult_SeqBAIJ_2; 2005 B->ops->multadd = MatMultAdd_SeqBAIJ_2; 2006 break; 2007 case 3: 2008 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 2009 B->ops->mult = MatMult_SeqBAIJ_3; 2010 B->ops->multadd = MatMultAdd_SeqBAIJ_3; 2011 break; 2012 case 4: 2013 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 2014 B->ops->mult = MatMult_SeqBAIJ_4; 2015 B->ops->multadd = MatMultAdd_SeqBAIJ_4; 2016 break; 2017 case 5: 2018 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 2019 B->ops->mult = MatMult_SeqBAIJ_5; 2020 B->ops->multadd = MatMultAdd_SeqBAIJ_5; 2021 break; 2022 case 6: 2023 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6; 2024 B->ops->mult = MatMult_SeqBAIJ_6; 2025 B->ops->multadd = MatMultAdd_SeqBAIJ_6; 2026 break; 2027 case 7: 2028 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7; 2029 B->ops->mult = MatMult_SeqBAIJ_7; 2030 B->ops->multadd = MatMultAdd_SeqBAIJ_7; 2031 break; 2032 default: 2033 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N; 2034 B->ops->mult = MatMult_SeqBAIJ_N; 2035 B->ops->multadd = MatMultAdd_SeqBAIJ_N; 2036 break; 2037 } 2038 } 2039 b->bs = bs; 2040 b->mbs = mbs; 2041 b->nbs = nbs; 2042 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2043 if (!nnz) { 2044 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2045 else if (nz <= 0) nz = 1; 2046 for (i=0; i<mbs; i++) b->imax[i] = nz; 2047 nz = nz*mbs; 2048 } else { 2049 nz = 0; 2050 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2051 } 2052 2053 /* allocate the matrix space */ 2054 len = nz*sizeof(int) + nz*bs2*sizeof(MatScalar) + (B->m+1)*sizeof(int); 2055 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2056 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2057 b->j = (int*)(b->a + nz*bs2); 2058 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2059 b->i = b->j + nz; 2060 b->singlemalloc = PETSC_TRUE; 2061 2062 b->i[0] = 0; 2063 for (i=1; i<mbs+1; i++) { 2064 b->i[i] = b->i[i-1] + b->imax[i-1]; 2065 } 2066 2067 /* b->ilen will count nonzeros in each block row so far. */ 2068 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2069 PetscLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 2070 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 2071 2072 b->bs = bs; 2073 b->bs2 = bs2; 2074 b->mbs = mbs; 2075 b->nz = 0; 2076 b->maxnz = nz*bs2; 2077 B->info.nz_unneeded = (PetscReal)b->maxnz; 2078 PetscFunctionReturn(0); 2079 } 2080 2081