1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: baij.c,v 1.156 1999/01/27 19:47:47 bsmith Exp curfman $"; 3 #endif 4 5 /* 6 Defines the basic matrix operations for the BAIJ (compressed row) 7 matrix storage format. 8 */ 9 #include "sys.h" 10 #include "src/mat/impls/baij/seq/baij.h" 11 #include "src/vec/vecimpl.h" 12 #include "src/inline/spops.h" 13 14 #define CHUNKSIZE 10 15 16 /* 17 Checks for missing diagonals 18 */ 19 #undef __FUNC__ 20 #define __FUNC__ "MatMissingDiag_SeqBAIJ" 21 int MatMissingDiag_SeqBAIJ(Mat A) 22 { 23 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 24 int *diag = a->diag, *jj = a->j,i; 25 26 PetscFunctionBegin; 27 for ( i=0; i<a->mbs; i++ ) { 28 if (jj[diag[i]] != i) { 29 SETERRQ1(1,1,"Matrix is missing diagonal number %d",i); 30 } 31 } 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNC__ 36 #define __FUNC__ "MatMarkDiag_SeqBAIJ" 37 int MatMarkDiag_SeqBAIJ(Mat A) 38 { 39 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 40 int i,j, *diag, m = a->mbs; 41 42 PetscFunctionBegin; 43 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 44 PLogObjectMemory(A,(m+1)*sizeof(int)); 45 for ( i=0; i<m; i++ ) { 46 diag[i] = a->i[i+1]; 47 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 48 if (a->j[j] == i) { 49 diag[i] = j; 50 break; 51 } 52 } 53 } 54 a->diag = diag; 55 PetscFunctionReturn(0); 56 } 57 58 59 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 60 61 #undef __FUNC__ 62 #define __FUNC__ "MatGetRowIJ_SeqBAIJ" 63 static int MatGetRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 64 PetscTruth *done) 65 { 66 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 67 int ierr, n = a->mbs,i; 68 69 PetscFunctionBegin; 70 *nn = n; 71 if (!ia) PetscFunctionReturn(0); 72 if (symmetric) { 73 ierr = MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 74 } else if (oshift == 1) { 75 /* temporarily add 1 to i and j indices */ 76 int nz = a->i[n] + 1; 77 for ( i=0; i<nz; i++ ) a->j[i]++; 78 for ( i=0; i<n+1; i++ ) a->i[i]++; 79 *ia = a->i; *ja = a->j; 80 } else { 81 *ia = a->i; *ja = a->j; 82 } 83 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNC__ 88 #define __FUNC__ "MatRestoreRowIJ_SeqBAIJ" 89 static int MatRestoreRowIJ_SeqBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 90 PetscTruth *done) 91 { 92 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 93 int i,n = a->mbs; 94 95 PetscFunctionBegin; 96 if (!ia) PetscFunctionReturn(0); 97 if (symmetric) { 98 PetscFree(*ia); 99 PetscFree(*ja); 100 } else if (oshift == 1) { 101 int nz = a->i[n]; 102 for ( i=0; i<nz; i++ ) a->j[i]--; 103 for ( i=0; i<n+1; i++ ) a->i[i]--; 104 } 105 PetscFunctionReturn(0); 106 } 107 108 #undef __FUNC__ 109 #define __FUNC__ "MatGetBlockSize_SeqBAIJ" 110 int MatGetBlockSize_SeqBAIJ(Mat mat, int *bs) 111 { 112 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *) mat->data; 113 114 PetscFunctionBegin; 115 *bs = baij->bs; 116 PetscFunctionReturn(0); 117 } 118 119 120 #undef __FUNC__ 121 #define __FUNC__ "MatDestroy_SeqBAIJ" 122 int MatDestroy_SeqBAIJ(Mat A) 123 { 124 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 125 int ierr; 126 127 if (--A->refct > 0) PetscFunctionReturn(0); 128 129 if (A->mapping) { 130 ierr = ISLocalToGlobalMappingDestroy(A->mapping); CHKERRQ(ierr); 131 } 132 if (A->bmapping) { 133 ierr = ISLocalToGlobalMappingDestroy(A->bmapping); CHKERRQ(ierr); 134 } 135 if (A->rmap) { 136 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 137 } 138 if (A->cmap) { 139 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 140 } 141 #if defined(USE_PETSC_LOG) 142 PLogObjectState((PetscObject) A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 143 #endif 144 PetscFree(a->a); 145 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 146 if (a->diag) PetscFree(a->diag); 147 if (a->ilen) PetscFree(a->ilen); 148 if (a->imax) PetscFree(a->imax); 149 if (a->solve_work) PetscFree(a->solve_work); 150 if (a->mult_work) PetscFree(a->mult_work); 151 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 152 PetscFree(a); 153 PLogObjectDestroy(A); 154 PetscHeaderDestroy(A); 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNC__ 159 #define __FUNC__ "MatSetOption_SeqBAIJ" 160 int MatSetOption_SeqBAIJ(Mat A,MatOption op) 161 { 162 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 163 164 PetscFunctionBegin; 165 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 166 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 167 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 168 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 169 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 170 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 171 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 172 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 173 else if (op == MAT_ROWS_SORTED || 174 op == MAT_ROWS_UNSORTED || 175 op == MAT_SYMMETRIC || 176 op == MAT_STRUCTURALLY_SYMMETRIC || 177 op == MAT_YES_NEW_DIAGONALS || 178 op == MAT_IGNORE_OFF_PROC_ENTRIES || 179 op == MAT_USE_HASH_TABLE) { 180 PLogInfo(A,"MatSetOption_SeqBAIJ:Option ignored\n"); 181 } else if (op == MAT_NO_NEW_DIAGONALS) { 182 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 183 } else { 184 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 185 } 186 PetscFunctionReturn(0); 187 } 188 189 190 #undef __FUNC__ 191 #define __FUNC__ "MatGetSize_SeqBAIJ" 192 int MatGetSize_SeqBAIJ(Mat A,int *m,int *n) 193 { 194 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 195 196 PetscFunctionBegin; 197 if (m) *m = a->m; 198 if (n) *n = a->n; 199 PetscFunctionReturn(0); 200 } 201 202 #undef __FUNC__ 203 #define __FUNC__ "MatGetOwnershipRange_SeqBAIJ" 204 int MatGetOwnershipRange_SeqBAIJ(Mat A,int *m,int *n) 205 { 206 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 207 208 PetscFunctionBegin; 209 *m = 0; *n = a->m; 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNC__ 214 #define __FUNC__ "MatGetRow_SeqBAIJ" 215 int MatGetRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 216 { 217 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 218 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*idx_i,bs2; 219 MatScalar *aa,*aa_i; 220 Scalar *v_i; 221 222 PetscFunctionBegin; 223 bs = a->bs; 224 ai = a->i; 225 aj = a->j; 226 aa = a->a; 227 bs2 = a->bs2; 228 229 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 230 231 bn = row/bs; /* Block number */ 232 bp = row % bs; /* Block Position */ 233 M = ai[bn+1] - ai[bn]; 234 *nz = bs*M; 235 236 if (v) { 237 *v = 0; 238 if (*nz) { 239 *v = (Scalar *) PetscMalloc( (*nz)*sizeof(Scalar) ); CHKPTRQ(*v); 240 for ( i=0; i<M; i++ ) { /* for each block in the block row */ 241 v_i = *v + i*bs; 242 aa_i = aa + bs2*(ai[bn] + i); 243 for ( j=bp,k=0; j<bs2; j+=bs,k++ ) {v_i[k] = aa_i[j];} 244 } 245 } 246 } 247 248 if (idx) { 249 *idx = 0; 250 if (*nz) { 251 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 252 for ( i=0; i<M; i++ ) { /* for each block in the block row */ 253 idx_i = *idx + i*bs; 254 itmp = bs*aj[ai[bn] + i]; 255 for ( j=0; j<bs; j++ ) {idx_i[j] = itmp++;} 256 } 257 } 258 } 259 PetscFunctionReturn(0); 260 } 261 262 #undef __FUNC__ 263 #define __FUNC__ "MatRestoreRow_SeqBAIJ" 264 int MatRestoreRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 265 { 266 PetscFunctionBegin; 267 if (idx) {if (*idx) PetscFree(*idx);} 268 if (v) {if (*v) PetscFree(*v);} 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNC__ 273 #define __FUNC__ "MatTranspose_SeqBAIJ" 274 int MatTranspose_SeqBAIJ(Mat A,Mat *B) 275 { 276 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 277 Mat C; 278 int i,j,k,ierr,*aj=a->j,*ai=a->i,bs=a->bs,mbs=a->mbs,nbs=a->nbs,len,*col; 279 int *rows,*cols,bs2=a->bs2; 280 MatScalar *array = a->a; 281 282 PetscFunctionBegin; 283 if (B==PETSC_NULL && mbs!=nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Square matrix only for in-place"); 284 col = (int *) PetscMalloc((1+nbs)*sizeof(int)); CHKPTRQ(col); 285 PetscMemzero(col,(1+nbs)*sizeof(int)); 286 287 for ( i=0; i<ai[mbs]; i++ ) col[aj[i]] += 1; 288 ierr = MatCreateSeqBAIJ(A->comm,bs,a->n,a->m,PETSC_NULL,col,&C); CHKERRQ(ierr); 289 PetscFree(col); 290 rows = (int *) PetscMalloc(2*bs*sizeof(int)); CHKPTRQ(rows); 291 cols = rows + bs; 292 for ( i=0; i<mbs; i++ ) { 293 cols[0] = i*bs; 294 for (k=1; k<bs; k++ ) cols[k] = cols[k-1] + 1; 295 len = ai[i+1] - ai[i]; 296 for ( j=0; j<len; j++ ) { 297 rows[0] = (*aj++)*bs; 298 for (k=1; k<bs; k++ ) rows[k] = rows[k-1] + 1; 299 ierr = MatSetValues(C,bs,rows,bs,cols,array,INSERT_VALUES); CHKERRQ(ierr); 300 array += bs2; 301 } 302 } 303 PetscFree(rows); 304 305 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 306 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 307 308 if (B != PETSC_NULL) { 309 *B = C; 310 } else { 311 PetscOps *Abops; 312 MatOps Aops; 313 314 /* This isn't really an in-place transpose */ 315 PetscFree(a->a); 316 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 317 if (a->diag) PetscFree(a->diag); 318 if (a->ilen) PetscFree(a->ilen); 319 if (a->imax) PetscFree(a->imax); 320 if (a->solve_work) PetscFree(a->solve_work); 321 PetscFree(a); 322 323 324 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 325 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 326 327 /* 328 This is horrible, horrible code. We need to keep the 329 A pointers for the bops and ops but copy everything 330 else from C. 331 */ 332 Abops = A->bops; 333 Aops = A->ops; 334 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 335 A->bops = Abops; 336 A->ops = Aops; 337 A->qlist = 0; 338 339 PetscHeaderDestroy(C); 340 } 341 PetscFunctionReturn(0); 342 } 343 344 345 346 347 #undef __FUNC__ 348 #define __FUNC__ "MatView_SeqBAIJ_Binary" 349 static int MatView_SeqBAIJ_Binary(Mat A,Viewer viewer) 350 { 351 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 352 int i, fd, *col_lens, ierr, bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 353 Scalar *aa; 354 355 PetscFunctionBegin; 356 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 357 col_lens = (int *) PetscMalloc((4+a->m)*sizeof(int));CHKPTRQ(col_lens); 358 col_lens[0] = MAT_COOKIE; 359 360 col_lens[1] = a->m; 361 col_lens[2] = a->n; 362 col_lens[3] = a->nz*bs2; 363 364 /* store lengths of each row and write (including header) to file */ 365 count = 0; 366 for ( i=0; i<a->mbs; i++ ) { 367 for ( j=0; j<bs; j++ ) { 368 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 369 } 370 } 371 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1); CHKERRQ(ierr); 372 PetscFree(col_lens); 373 374 /* store column indices (zero start index) */ 375 jj = (int *) PetscMalloc( (a->nz+1)*bs2*sizeof(int) ); CHKPTRQ(jj); 376 count = 0; 377 for ( i=0; i<a->mbs; i++ ) { 378 for ( j=0; j<bs; j++ ) { 379 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 380 for ( l=0; l<bs; l++ ) { 381 jj[count++] = bs*a->j[k] + l; 382 } 383 } 384 } 385 } 386 ierr = PetscBinaryWrite(fd,jj,bs2*a->nz,PETSC_INT,0); CHKERRQ(ierr); 387 PetscFree(jj); 388 389 /* store nonzero values */ 390 aa = (Scalar *) PetscMalloc((a->nz+1)*bs2*sizeof(Scalar)); CHKPTRQ(aa); 391 count = 0; 392 for ( i=0; i<a->mbs; i++ ) { 393 for ( j=0; j<bs; j++ ) { 394 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 395 for ( l=0; l<bs; l++ ) { 396 aa[count++] = a->a[bs2*k + l*bs + j]; 397 } 398 } 399 } 400 } 401 ierr = PetscBinaryWrite(fd,aa,bs2*a->nz,PETSC_SCALAR,0); CHKERRQ(ierr); 402 PetscFree(aa); 403 PetscFunctionReturn(0); 404 } 405 406 #undef __FUNC__ 407 #define __FUNC__ "MatView_SeqBAIJ_ASCII" 408 static int MatView_SeqBAIJ_ASCII(Mat A,Viewer viewer) 409 { 410 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 411 int ierr, i,j,format,bs = a->bs,k,l,bs2=a->bs2; 412 FILE *fd; 413 char *outputname; 414 415 PetscFunctionBegin; 416 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 417 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 418 ierr = ViewerGetFormat(viewer,&format); 419 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 420 ViewerASCIIPrintf(viewer," block size is %d\n",bs); 421 } else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 422 SETERRQ(PETSC_ERR_SUP,0,"Socket format not supported"); 423 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 424 for ( i=0; i<a->mbs; i++ ) { 425 for ( j=0; j<bs; j++ ) { 426 fprintf(fd,"row %d:",i*bs+j); 427 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 428 for ( l=0; l<bs; l++ ) { 429 #if defined(USE_PETSC_COMPLEX) 430 if (PetscImaginary(a->a[bs2*k + l*bs + j]) > 0.0 && PetscReal(a->a[bs2*k + l*bs + j]) != 0.0) { 431 fprintf(fd," %d %g + %g i",bs*a->j[k]+l, 432 PetscReal(a->a[bs2*k + l*bs + j]),PetscImaginary(a->a[bs2*k + l*bs + j])); 433 } else if (PetscImaginary(a->a[bs2*k + l*bs + j]) < 0.0 && PetscReal(a->a[bs2*k + l*bs + j]) != 0.0) { 434 fprintf(fd," %d %g - %g i",bs*a->j[k]+l, 435 PetscReal(a->a[bs2*k + l*bs + j]),-PetscImaginary(a->a[bs2*k + l*bs + j])); 436 } else if (PetscReal(a->a[bs2*k + l*bs + j]) != 0.0) { 437 fprintf(fd," %d %g ",bs*a->j[k]+l,PetscReal(a->a[bs2*k + l*bs + j])); 438 } 439 #else 440 if (a->a[bs2*k + l*bs + j] != 0.0) { 441 fprintf(fd," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]); 442 } 443 #endif 444 } 445 } 446 fprintf(fd,"\n"); 447 } 448 } 449 } else { 450 for ( i=0; i<a->mbs; i++ ) { 451 for ( j=0; j<bs; j++ ) { 452 fprintf(fd,"row %d:",i*bs+j); 453 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 454 for ( l=0; l<bs; l++ ) { 455 #if defined(USE_PETSC_COMPLEX) 456 if (PetscImaginary(a->a[bs2*k + l*bs + j]) > 0.0) { 457 fprintf(fd," %d %g + %g i",bs*a->j[k]+l, 458 PetscReal(a->a[bs2*k + l*bs + j]),PetscImaginary(a->a[bs2*k + l*bs + j])); 459 } else if (PetscImaginary(a->a[bs2*k + l*bs + j]) < 0.0) { 460 fprintf(fd," %d %g - %g i",bs*a->j[k]+l, 461 PetscReal(a->a[bs2*k + l*bs + j]),-PetscImaginary(a->a[bs2*k + l*bs + j])); 462 } else { 463 fprintf(fd," %d %g ",bs*a->j[k]+l,PetscReal(a->a[bs2*k + l*bs + j])); 464 } 465 #else 466 fprintf(fd," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]); 467 #endif 468 } 469 } 470 fprintf(fd,"\n"); 471 } 472 } 473 } 474 fflush(fd); 475 PetscFunctionReturn(0); 476 } 477 478 #undef __FUNC__ 479 #define __FUNC__ "MatView_SeqBAIJ_Draw_Zoom" 480 static int MatView_SeqBAIJ_Draw_Zoom(Draw draw,void *Aa) 481 { 482 Mat A = (Mat) Aa; 483 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *) A->data; 484 int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2,rank; 485 double xl,yl,xr,yr,x_l,x_r,y_l,y_r; 486 MatScalar *aa; 487 MPI_Comm comm; 488 Viewer viewer; 489 490 PetscFunctionBegin; 491 /* 492 This is nasty. If this is called from an originally parallel matrix 493 then all processes call this, but only the first has the matrix so the 494 rest should return immediately. 495 */ 496 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 497 MPI_Comm_rank(comm,&rank); 498 if (rank) PetscFunctionReturn(0); 499 500 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*) &viewer);CHKERRQ(ierr); 501 502 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr); CHKERRQ(ierr); 503 504 /* loop over matrix elements drawing boxes */ 505 color = DRAW_BLUE; 506 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 507 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 508 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 509 x_l = a->j[j]*bs; x_r = x_l + 1.0; 510 aa = a->a + j*bs2; 511 for ( k=0; k<bs; k++ ) { 512 for ( l=0; l<bs; l++ ) { 513 if (PetscReal(*aa++) >= 0.) continue; 514 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 515 } 516 } 517 } 518 } 519 color = DRAW_CYAN; 520 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 521 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 522 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 523 x_l = a->j[j]*bs; x_r = x_l + 1.0; 524 aa = a->a + j*bs2; 525 for ( k=0; k<bs; k++ ) { 526 for ( l=0; l<bs; l++ ) { 527 if (PetscReal(*aa++) != 0.) continue; 528 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 529 } 530 } 531 } 532 } 533 534 color = DRAW_RED; 535 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 536 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 537 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 538 x_l = a->j[j]*bs; x_r = x_l + 1.0; 539 aa = a->a + j*bs2; 540 for ( k=0; k<bs; k++ ) { 541 for ( l=0; l<bs; l++ ) { 542 if (PetscReal(*aa++) <= 0.) continue; 543 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 544 } 545 } 546 } 547 } 548 PetscFunctionReturn(0); 549 } 550 551 #undef __FUNC__ 552 #define __FUNC__ "MatView_SeqBAIJ_Draw" 553 static int MatView_SeqBAIJ_Draw(Mat A,Viewer viewer) 554 { 555 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *) A->data; 556 int ierr; 557 double xl,yl,xr,yr,w,h; 558 Draw draw; 559 PetscTruth isnull; 560 561 PetscFunctionBegin; 562 563 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 564 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 565 566 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 567 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 568 xr += w; yr += h; xl = -w; yl = -h; 569 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 570 ierr = DrawZoom(draw,MatView_SeqBAIJ_Draw_Zoom,A); CHKERRQ(ierr); 571 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 572 PetscFunctionReturn(0); 573 } 574 575 #undef __FUNC__ 576 #define __FUNC__ "MatView_SeqBAIJ" 577 int MatView_SeqBAIJ(Mat A,Viewer viewer) 578 { 579 ViewerType vtype; 580 int ierr; 581 582 PetscFunctionBegin; 583 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 584 if (PetscTypeCompare(vtype,SOCKET_VIEWER)) { 585 SETERRQ(PETSC_ERR_SUP,0,"Socket viewer not supported"); 586 } else if (PetscTypeCompare(vtype,ASCII_VIEWER)){ 587 ierr = MatView_SeqBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 588 } else if (PetscTypeCompare(vtype,BINARY_VIEWER)) { 589 ierr = MatView_SeqBAIJ_Binary(A,viewer);CHKERRQ(ierr); 590 } else if (PetscTypeCompare(vtype,DRAW_VIEWER)) { 591 ierr = MatView_SeqBAIJ_Draw(A,viewer);CHKERRQ(ierr); 592 } else { 593 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 594 } 595 PetscFunctionReturn(0); 596 } 597 598 599 #undef __FUNC__ 600 #define __FUNC__ "MatGetValues_SeqBAIJ" 601 int MatGetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 602 { 603 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 604 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 605 int *ai = a->i, *ailen = a->ilen; 606 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 607 MatScalar *ap, *aa = a->a, zero = 0.0; 608 609 PetscFunctionBegin; 610 for ( k=0; k<m; k++ ) { /* loop over rows */ 611 row = im[k]; brow = row/bs; 612 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 613 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 614 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 615 nrow = ailen[brow]; 616 for ( l=0; l<n; l++ ) { /* loop over columns */ 617 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 618 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 619 col = in[l] ; 620 bcol = col/bs; 621 cidx = col%bs; 622 ridx = row%bs; 623 high = nrow; 624 low = 0; /* assume unsorted */ 625 while (high-low > 5) { 626 t = (low+high)/2; 627 if (rp[t] > bcol) high = t; 628 else low = t; 629 } 630 for ( i=low; i<high; i++ ) { 631 if (rp[i] > bcol) break; 632 if (rp[i] == bcol) { 633 *v++ = ap[bs2*i+bs*cidx+ridx]; 634 goto finished; 635 } 636 } 637 *v++ = zero; 638 finished:; 639 } 640 } 641 PetscFunctionReturn(0); 642 } 643 644 645 #undef __FUNC__ 646 #define __FUNC__ "MatSetValuesBlocked_SeqBAIJ" 647 int MatSetValuesBlocked_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 648 { 649 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 650 int *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 651 int *imax=a->imax,*ai=a->i,*ailen=a->ilen, roworiented=a->roworiented; 652 int *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=a->bs,stepval; 653 Scalar *value = v; 654 MatScalar *ap,*aa=a->a,*bap; 655 656 PetscFunctionBegin; 657 if (roworiented) { 658 stepval = (n-1)*bs; 659 } else { 660 stepval = (m-1)*bs; 661 } 662 for ( k=0; k<m; k++ ) { /* loop over added rows */ 663 row = im[k]; 664 if (row < 0) continue; 665 #if defined(USE_PETSC_BOPT_g) 666 if (row >= a->mbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 667 #endif 668 rp = aj + ai[row]; 669 ap = aa + bs2*ai[row]; 670 rmax = imax[row]; 671 nrow = ailen[row]; 672 low = 0; 673 for ( l=0; l<n; l++ ) { /* loop over added columns */ 674 if (in[l] < 0) continue; 675 #if defined(USE_PETSC_BOPT_g) 676 if (in[l] >= a->nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 677 #endif 678 col = in[l]; 679 if (roworiented) { 680 value = v + k*(stepval+bs)*bs + l*bs; 681 } else { 682 value = v + l*(stepval+bs)*bs + k*bs; 683 } 684 if (!sorted) low = 0; high = nrow; 685 while (high-low > 7) { 686 t = (low+high)/2; 687 if (rp[t] > col) high = t; 688 else low = t; 689 } 690 for ( i=low; i<high; i++ ) { 691 if (rp[i] > col) break; 692 if (rp[i] == col) { 693 bap = ap + bs2*i; 694 if (roworiented) { 695 if (is == ADD_VALUES) { 696 for ( ii=0; ii<bs; ii++,value+=stepval ) { 697 for (jj=ii; jj<bs2; jj+=bs ) { 698 bap[jj] += *value++; 699 } 700 } 701 } else { 702 for ( ii=0; ii<bs; ii++,value+=stepval ) { 703 for (jj=ii; jj<bs2; jj+=bs ) { 704 bap[jj] = *value++; 705 } 706 } 707 } 708 } else { 709 if (is == ADD_VALUES) { 710 for ( ii=0; ii<bs; ii++,value+=stepval ) { 711 for (jj=0; jj<bs; jj++ ) { 712 *bap++ += *value++; 713 } 714 } 715 } else { 716 for ( ii=0; ii<bs; ii++,value+=stepval ) { 717 for (jj=0; jj<bs; jj++ ) { 718 *bap++ = *value++; 719 } 720 } 721 } 722 } 723 goto noinsert2; 724 } 725 } 726 if (nonew == 1) goto noinsert2; 727 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 728 if (nrow >= rmax) { 729 /* there is no extra room in row, therefore enlarge */ 730 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 731 MatScalar *new_a; 732 733 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 734 735 /* malloc new storage space */ 736 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 737 new_a = (MatScalar *) PetscMalloc( len ); CHKPTRQ(new_a); 738 new_j = (int *) (new_a + bs2*new_nz); 739 new_i = new_j + new_nz; 740 741 /* copy over old data into new slots */ 742 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 743 for ( ii=row+1; ii<a->mbs+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 744 PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int)); 745 len = (new_nz - CHUNKSIZE - ai[row] - nrow); 746 PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int)); 747 PetscMemcpy(new_a,aa,(ai[row]+nrow)*bs2*sizeof(MatScalar)); 748 PetscMemzero(new_a+bs2*(ai[row]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar)); 749 PetscMemcpy(new_a+bs2*(ai[row]+nrow+CHUNKSIZE),aa+bs2*(ai[row]+nrow),bs2*len*sizeof(MatScalar)); 750 /* free up old matrix storage */ 751 PetscFree(a->a); 752 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 753 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 754 a->singlemalloc = 1; 755 756 rp = aj + ai[row]; ap = aa + bs2*ai[row]; 757 rmax = imax[row] = imax[row] + CHUNKSIZE; 758 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 759 a->maxnz += bs2*CHUNKSIZE; 760 a->reallocs++; 761 a->nz++; 762 } 763 N = nrow++ - 1; 764 /* shift up all the later entries in this row */ 765 for ( ii=N; ii>=i; ii-- ) { 766 rp[ii+1] = rp[ii]; 767 PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar)); 768 } 769 if (N >= i) PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar)); 770 rp[i] = col; 771 bap = ap + bs2*i; 772 if (roworiented) { 773 for ( ii=0; ii<bs; ii++,value+=stepval) { 774 for (jj=ii; jj<bs2; jj+=bs ) { 775 bap[jj] = *value++; 776 } 777 } 778 } else { 779 for ( ii=0; ii<bs; ii++,value+=stepval ) { 780 for (jj=0; jj<bs; jj++ ) { 781 *bap++ = *value++; 782 } 783 } 784 } 785 noinsert2:; 786 low = i; 787 } 788 ailen[row] = nrow; 789 } 790 PetscFunctionReturn(0); 791 } 792 793 794 #undef __FUNC__ 795 #define __FUNC__ "MatAssemblyEnd_SeqBAIJ" 796 int MatAssemblyEnd_SeqBAIJ(Mat A,MatAssemblyType mode) 797 { 798 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 799 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax; 800 int m = a->m,*ip, N, *ailen = a->ilen; 801 int mbs = a->mbs, bs2 = a->bs2,rmax = 0; 802 MatScalar *aa = a->a, *ap; 803 804 PetscFunctionBegin; 805 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 806 807 if (m) rmax = ailen[0]; 808 for ( i=1; i<mbs; i++ ) { 809 /* move each row back by the amount of empty slots (fshift) before it*/ 810 fshift += imax[i-1] - ailen[i-1]; 811 rmax = PetscMax(rmax,ailen[i]); 812 if (fshift) { 813 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 814 N = ailen[i]; 815 for ( j=0; j<N; j++ ) { 816 ip[j-fshift] = ip[j]; 817 PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar)); 818 } 819 } 820 ai[i] = ai[i-1] + ailen[i-1]; 821 } 822 if (mbs) { 823 fshift += imax[mbs-1] - ailen[mbs-1]; 824 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 825 } 826 /* reset ilen and imax for each row */ 827 for ( i=0; i<mbs; i++ ) { 828 ailen[i] = imax[i] = ai[i+1] - ai[i]; 829 } 830 a->nz = ai[mbs]; 831 832 /* diagonals may have moved, so kill the diagonal pointers */ 833 if (fshift && a->diag) { 834 PetscFree(a->diag); 835 PLogObjectMemory(A,-(m+1)*sizeof(int)); 836 a->diag = 0; 837 } 838 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n", 839 m,a->n,a->bs,fshift*bs2,a->nz*bs2); 840 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Number of mallocs during MatSetValues is %d\n", 841 a->reallocs); 842 PLogInfo(A,"MatAssemblyEnd_SeqBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 843 a->reallocs = 0; 844 A->info.nz_unneeded = (double)fshift*bs2; 845 846 PetscFunctionReturn(0); 847 } 848 849 850 /* idx should be of length atlease bs */ 851 #undef __FUNC__ 852 #define __FUNC__ "MatZeroRows_SeqBAIJ_Check_Block" 853 static int MatZeroRows_SeqBAIJ_Check_Block(int *idx, int bs, PetscTruth *flg) 854 { 855 int i,row; 856 857 PetscFunctionBegin; 858 row = idx[0]; 859 if (row%bs!=0) { *flg = PETSC_FALSE; PetscFunctionReturn(0); } 860 861 for ( i=1; i<bs; i++ ) { 862 if (row+i != idx[i]) { *flg = PETSC_FALSE; PetscFunctionReturn(0); } 863 } 864 *flg = PETSC_TRUE; 865 PetscFunctionReturn(0); 866 } 867 868 #undef __FUNC__ 869 #define __FUNC__ "MatZeroRows_SeqBAIJ" 870 int MatZeroRows_SeqBAIJ(Mat A,IS is, Scalar *diag) 871 { 872 Mat_SeqBAIJ *baij=(Mat_SeqBAIJ*)A->data; 873 IS is_local; 874 int ierr,i,j,count,m=baij->m,is_n,*is_idx,*rows,bs=baij->bs,bs2=baij->bs2; 875 PetscTruth flg; 876 Scalar zero = 0.0; 877 MatScalar *aa; 878 879 PetscFunctionBegin; 880 /* Make a copy of the IS and sort it */ 881 ierr = ISGetSize(is,&is_n);CHKERRQ(ierr); 882 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 883 ierr = ISCreateGeneral(A->comm,is_n,is_idx,&is_local); CHKERRQ(ierr); 884 ierr = ISSort(is_local); CHKERRQ(ierr); 885 ierr = ISGetIndices(is_local,&rows); CHKERRQ(ierr); 886 887 i = 0; 888 while (i < is_n) { 889 if (rows[i]<0 || rows[i]>m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 890 flg = PETSC_FALSE; 891 if (i+bs <= is_n) {ierr = MatZeroRows_SeqBAIJ_Check_Block(rows+i,bs,&flg); CHKERRQ(ierr); } 892 count = (baij->i[rows[i]/bs +1] - baij->i[rows[i]/bs])*bs; 893 aa = baij->a + baij->i[rows[i]/bs]*bs2 + (rows[i]%bs); 894 if (flg) { /* There exists a block of rows to be Zerowed */ 895 if (baij->ilen[rows[i]/bs] > 0) { 896 PetscMemzero(aa,count*bs*sizeof(MatScalar)); 897 baij->ilen[rows[i]/bs] = 1; 898 baij->j[baij->i[rows[i]/bs]] = rows[i]/bs; 899 } 900 i += bs; 901 } else { /* Zero out only the requested row */ 902 for ( j=0; j<count; j++ ) { 903 aa[0] = zero; 904 aa+=bs; 905 } 906 i++; 907 } 908 } 909 if (diag) { 910 for ( j=0; j<is_n; j++ ) { 911 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 912 } 913 } 914 ierr = ISRestoreIndices(is,&is_idx); CHKERRQ(ierr); 915 ierr = ISRestoreIndices(is_local,&rows); CHKERRQ(ierr); 916 ierr = ISDestroy(is_local); CHKERRQ(ierr); 917 ierr = MatAssemblyEnd_SeqBAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 918 PetscFunctionReturn(0); 919 } 920 921 #undef __FUNC__ 922 #define __FUNC__ "MatSetValues_SeqBAIJ" 923 int MatSetValues_SeqBAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 924 { 925 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 926 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 927 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 928 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 929 int ridx,cidx,bs2=a->bs2; 930 MatScalar *ap,value,*aa=a->a,*bap; 931 932 PetscFunctionBegin; 933 for ( k=0; k<m; k++ ) { /* loop over added rows */ 934 row = im[k]; brow = row/bs; 935 if (row < 0) continue; 936 #if defined(USE_PETSC_BOPT_g) 937 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 938 #endif 939 rp = aj + ai[brow]; 940 ap = aa + bs2*ai[brow]; 941 rmax = imax[brow]; 942 nrow = ailen[brow]; 943 low = 0; 944 for ( l=0; l<n; l++ ) { /* loop over added columns */ 945 if (in[l] < 0) continue; 946 #if defined(USE_PETSC_BOPT_g) 947 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 948 #endif 949 col = in[l]; bcol = col/bs; 950 ridx = row % bs; cidx = col % bs; 951 if (roworiented) { 952 value = v[l + k*n]; 953 } else { 954 value = v[k + l*m]; 955 } 956 if (!sorted) low = 0; high = nrow; 957 while (high-low > 7) { 958 t = (low+high)/2; 959 if (rp[t] > bcol) high = t; 960 else low = t; 961 } 962 for ( i=low; i<high; i++ ) { 963 if (rp[i] > bcol) break; 964 if (rp[i] == bcol) { 965 bap = ap + bs2*i + bs*cidx + ridx; 966 if (is == ADD_VALUES) *bap += value; 967 else *bap = value; 968 goto noinsert1; 969 } 970 } 971 if (nonew == 1) goto noinsert1; 972 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 973 if (nrow >= rmax) { 974 /* there is no extra room in row, therefore enlarge */ 975 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 976 MatScalar *new_a; 977 978 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 979 980 /* Malloc new storage space */ 981 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 982 new_a = (MatScalar *) PetscMalloc( len ); CHKPTRQ(new_a); 983 new_j = (int *) (new_a + bs2*new_nz); 984 new_i = new_j + new_nz; 985 986 /* copy over old data into new slots */ 987 for ( ii=0; ii<brow+1; ii++ ) {new_i[ii] = ai[ii];} 988 for ( ii=brow+1; ii<a->mbs+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 989 PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int)); 990 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 991 PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow, 992 len*sizeof(int)); 993 PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(MatScalar)); 994 PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar)); 995 PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE), 996 aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(MatScalar)); 997 /* free up old matrix storage */ 998 PetscFree(a->a); 999 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 1000 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 1001 a->singlemalloc = 1; 1002 1003 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 1004 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 1005 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 1006 a->maxnz += bs2*CHUNKSIZE; 1007 a->reallocs++; 1008 a->nz++; 1009 } 1010 N = nrow++ - 1; 1011 /* shift up all the later entries in this row */ 1012 for ( ii=N; ii>=i; ii-- ) { 1013 rp[ii+1] = rp[ii]; 1014 PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar)); 1015 } 1016 if (N>=i) PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar)); 1017 rp[i] = bcol; 1018 ap[bs2*i + bs*cidx + ridx] = value; 1019 noinsert1:; 1020 low = i; 1021 } 1022 ailen[brow] = nrow; 1023 } 1024 PetscFunctionReturn(0); 1025 } 1026 1027 extern int MatLUFactorSymbolic_SeqBAIJ(Mat,IS,IS,double,Mat*); 1028 extern int MatLUFactor_SeqBAIJ(Mat,IS,IS,double); 1029 extern int MatIncreaseOverlap_SeqBAIJ(Mat,int,IS*,int); 1030 extern int MatGetSubMatrix_SeqBAIJ(Mat,IS,IS,int,MatReuse,Mat*); 1031 extern int MatGetSubMatrices_SeqBAIJ(Mat,int,IS*,IS*,MatReuse,Mat**); 1032 extern int MatMultTrans_SeqBAIJ(Mat,Vec,Vec); 1033 extern int MatMultTransAdd_SeqBAIJ(Mat,Vec,Vec,Vec); 1034 extern int MatScale_SeqBAIJ(Scalar*,Mat); 1035 extern int MatNorm_SeqBAIJ(Mat,NormType,double *); 1036 extern int MatEqual_SeqBAIJ(Mat,Mat, PetscTruth*); 1037 extern int MatGetDiagonal_SeqBAIJ(Mat,Vec); 1038 extern int MatDiagonalScale_SeqBAIJ(Mat,Vec,Vec); 1039 extern int MatGetInfo_SeqBAIJ(Mat,MatInfoType,MatInfo *); 1040 extern int MatZeroEntries_SeqBAIJ(Mat); 1041 1042 extern int MatSolve_SeqBAIJ_N(Mat,Vec,Vec); 1043 extern int MatSolve_SeqBAIJ_1(Mat,Vec,Vec); 1044 extern int MatSolve_SeqBAIJ_2(Mat,Vec,Vec); 1045 extern int MatSolve_SeqBAIJ_3(Mat,Vec,Vec); 1046 extern int MatSolve_SeqBAIJ_4(Mat,Vec,Vec); 1047 extern int MatSolve_SeqBAIJ_5(Mat,Vec,Vec); 1048 extern int MatSolve_SeqBAIJ_7(Mat,Vec,Vec); 1049 1050 extern int MatLUFactorNumeric_SeqBAIJ_N(Mat,Mat*); 1051 extern int MatLUFactorNumeric_SeqBAIJ_1(Mat,Mat*); 1052 extern int MatLUFactorNumeric_SeqBAIJ_2(Mat,Mat*); 1053 extern int MatLUFactorNumeric_SeqBAIJ_3(Mat,Mat*); 1054 extern int MatLUFactorNumeric_SeqBAIJ_4(Mat,Mat*); 1055 extern int MatLUFactorNumeric_SeqBAIJ_5(Mat,Mat*); 1056 1057 extern int MatMult_SeqBAIJ_1(Mat,Vec,Vec); 1058 extern int MatMult_SeqBAIJ_2(Mat,Vec,Vec); 1059 extern int MatMult_SeqBAIJ_3(Mat,Vec,Vec); 1060 extern int MatMult_SeqBAIJ_4(Mat,Vec,Vec); 1061 extern int MatMult_SeqBAIJ_5(Mat,Vec,Vec); 1062 extern int MatMult_SeqBAIJ_7(Mat,Vec,Vec); 1063 extern int MatMult_SeqBAIJ_N(Mat,Vec,Vec); 1064 1065 extern int MatMultAdd_SeqBAIJ_1(Mat,Vec,Vec,Vec); 1066 extern int MatMultAdd_SeqBAIJ_2(Mat,Vec,Vec,Vec); 1067 extern int MatMultAdd_SeqBAIJ_3(Mat,Vec,Vec,Vec); 1068 extern int MatMultAdd_SeqBAIJ_4(Mat,Vec,Vec,Vec); 1069 extern int MatMultAdd_SeqBAIJ_5(Mat,Vec,Vec,Vec); 1070 extern int MatMultAdd_SeqBAIJ_7(Mat,Vec,Vec,Vec); 1071 extern int MatMultAdd_SeqBAIJ_N(Mat,Vec,Vec,Vec); 1072 1073 #undef __FUNC__ 1074 #define __FUNC__ "MatILUFactor_SeqBAIJ" 1075 int MatILUFactor_SeqBAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1076 { 1077 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) inA->data; 1078 Mat outA; 1079 int ierr; 1080 PetscTruth row_identity, col_identity; 1081 1082 PetscFunctionBegin; 1083 if (info && info->fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1084 ierr = ISIdentity(row,&row_identity); CHKERRQ(ierr); 1085 ierr = ISIdentity(col,&col_identity); CHKERRQ(ierr); 1086 if (!row_identity || !col_identity) { 1087 SETERRQ(1,1,"Row and column permutations must be identity"); 1088 } 1089 1090 outA = inA; 1091 inA->factor = FACTOR_LU; 1092 a->row = row; 1093 a->col = col; 1094 1095 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1096 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1097 PLogObjectParent(inA,a->icol); 1098 1099 if (!a->solve_work) { 1100 a->solve_work = (Scalar *) PetscMalloc((a->m+a->bs)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1101 PLogObjectMemory(inA,(a->m+a->bs)*sizeof(Scalar)); 1102 } 1103 1104 if (!a->diag) { 1105 ierr = MatMarkDiag_SeqBAIJ(inA); CHKERRQ(ierr); 1106 } 1107 /* 1108 Blocksize 4 and 5 have a special faster factorization/solver for ILU(0) factorization 1109 with natural ordering 1110 */ 1111 if (a->bs == 4) { 1112 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 1113 inA->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 1114 PLogInfo(inA,"Using special in-place natural ordering factor and solve BS=4\n"); 1115 } else if (a->bs == 5) { 1116 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 1117 inA->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering; 1118 PLogInfo(inA,"Using special in-place natural ordering factor and solve BS=5\n"); 1119 } 1120 1121 ierr = MatLUFactorNumeric(inA,&outA); CHKERRQ(ierr); 1122 1123 1124 PetscFunctionReturn(0); 1125 } 1126 #undef __FUNC__ 1127 #define __FUNC__ "MatPrintHelp_SeqBAIJ" 1128 int MatPrintHelp_SeqBAIJ(Mat A) 1129 { 1130 static int called = 0; 1131 MPI_Comm comm = A->comm; 1132 1133 PetscFunctionBegin; 1134 if (called) {PetscFunctionReturn(0);} else called = 1; 1135 (*PetscHelpPrintf)(comm," Options for MATSEQBAIJ and MATMPIBAIJ matrix formats (the defaults):\n"); 1136 (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n"); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 EXTERN_C_BEGIN 1141 #undef __FUNC__ 1142 #define __FUNC__ "MatSeqBAIJSetColumnIndices_SeqBAIJ" 1143 int MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat,int *indices) 1144 { 1145 Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data; 1146 int i,nz,n; 1147 1148 PetscFunctionBegin; 1149 nz = baij->maxnz; 1150 n = baij->n; 1151 for (i=0; i<nz; i++) { 1152 baij->j[i] = indices[i]; 1153 } 1154 baij->nz = nz; 1155 for ( i=0; i<n; i++ ) { 1156 baij->ilen[i] = baij->imax[i]; 1157 } 1158 1159 PetscFunctionReturn(0); 1160 } 1161 EXTERN_C_END 1162 1163 #undef __FUNC__ 1164 #define __FUNC__ "MatSeqBAIJSetColumnIndices" 1165 /*@ 1166 MatSeqBAIJSetColumnIndices - Set the column indices for all the rows 1167 in the matrix. 1168 1169 Input Parameters: 1170 + mat - the SeqBAIJ matrix 1171 - indices - the column indices 1172 1173 Notes: 1174 This can be called if you have precomputed the nonzero structure of the 1175 matrix and want to provide it to the matrix object to improve the performance 1176 of the MatSetValues() operation. 1177 1178 You MUST have set the correct numbers of nonzeros per row in the call to 1179 MatCreateSeqBAIJ(). 1180 1181 MUST be called before any calls to MatSetValues(); 1182 1183 @*/ 1184 int MatSeqBAIJSetColumnIndices(Mat mat,int *indices) 1185 { 1186 int ierr,(*f)(Mat,int *); 1187 1188 PetscFunctionBegin; 1189 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1190 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqBAIJSetColumnIndices_C",(void **)&f); 1191 CHKERRQ(ierr); 1192 if (f) { 1193 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1194 } else { 1195 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1196 } 1197 PetscFunctionReturn(0); 1198 } 1199 1200 /* -------------------------------------------------------------------*/ 1201 static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ, 1202 MatGetRow_SeqBAIJ, 1203 MatRestoreRow_SeqBAIJ, 1204 MatMult_SeqBAIJ_N, 1205 MatMultAdd_SeqBAIJ_N, 1206 MatMultTrans_SeqBAIJ, 1207 MatMultTransAdd_SeqBAIJ, 1208 MatSolve_SeqBAIJ_N, 1209 0, 1210 0, 1211 0, 1212 MatLUFactor_SeqBAIJ, 1213 0, 1214 0, 1215 MatTranspose_SeqBAIJ, 1216 MatGetInfo_SeqBAIJ, 1217 MatEqual_SeqBAIJ, 1218 MatGetDiagonal_SeqBAIJ, 1219 MatDiagonalScale_SeqBAIJ, 1220 MatNorm_SeqBAIJ, 1221 0, 1222 MatAssemblyEnd_SeqBAIJ, 1223 0, 1224 MatSetOption_SeqBAIJ, 1225 MatZeroEntries_SeqBAIJ, 1226 MatZeroRows_SeqBAIJ, 1227 MatLUFactorSymbolic_SeqBAIJ, 1228 MatLUFactorNumeric_SeqBAIJ_N, 1229 0, 1230 0, 1231 MatGetSize_SeqBAIJ, 1232 MatGetSize_SeqBAIJ, 1233 MatGetOwnershipRange_SeqBAIJ, 1234 MatILUFactorSymbolic_SeqBAIJ, 1235 0, 1236 0, 1237 0, 1238 MatDuplicate_SeqBAIJ, 1239 0, 1240 0, 1241 MatILUFactor_SeqBAIJ, 1242 0, 1243 0, 1244 MatGetSubMatrices_SeqBAIJ, 1245 MatIncreaseOverlap_SeqBAIJ, 1246 MatGetValues_SeqBAIJ, 1247 0, 1248 MatPrintHelp_SeqBAIJ, 1249 MatScale_SeqBAIJ, 1250 0, 1251 0, 1252 0, 1253 MatGetBlockSize_SeqBAIJ, 1254 MatGetRowIJ_SeqBAIJ, 1255 MatRestoreRowIJ_SeqBAIJ, 1256 0, 1257 0, 1258 0, 1259 0, 1260 0, 1261 0, 1262 MatSetValuesBlocked_SeqBAIJ, 1263 MatGetSubMatrix_SeqBAIJ, 1264 0, 1265 0, 1266 MatGetMaps_Petsc}; 1267 1268 #undef __FUNC__ 1269 #define __FUNC__ "MatCreateSeqBAIJ" 1270 /*@C 1271 MatCreateSeqBAIJ - Creates a sparse matrix in block AIJ (block 1272 compressed row) format. For good matrix assembly performance the 1273 user should preallocate the matrix storage by setting the parameter nz 1274 (or the array nzz). By setting these parameters accurately, performance 1275 during matrix assembly can be increased by more than a factor of 50. 1276 1277 Collective on MPI_Comm 1278 1279 Input Parameters: 1280 + comm - MPI communicator, set to PETSC_COMM_SELF 1281 . bs - size of block 1282 . m - number of rows 1283 . n - number of columns 1284 . nz - number of block nonzeros per block row (same for all rows) 1285 - nzz - array containing the number of block nonzeros in the various block rows 1286 (possibly different for each block row) or PETSC_NULL 1287 1288 Output Parameter: 1289 . A - the matrix 1290 1291 Options Database Keys: 1292 . -mat_no_unroll - uses code that does not unroll the loops in the 1293 block calculations (much slower) 1294 . -mat_block_size - size of the blocks to use 1295 1296 Notes: 1297 The block AIJ format is fully compatible with standard Fortran 77 1298 storage. That is, the stored row and column indices can begin at 1299 either one (as in Fortran) or zero. See the users' manual for details. 1300 1301 Specify the preallocated storage with either nz or nnz (not both). 1302 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1303 allocation. For additional details, see the users manual chapter on 1304 matrices. 1305 1306 Level: intermediate 1307 1308 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPIBAIJ() 1309 @*/ 1310 int MatCreateSeqBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz, Mat *A) 1311 { 1312 Mat B; 1313 Mat_SeqBAIJ *b; 1314 int i,len,ierr,flg,mbs=m/bs,nbs=n/bs,bs2=bs*bs,size; 1315 1316 PetscFunctionBegin; 1317 MPI_Comm_size(comm,&size); 1318 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1319 1320 ierr = OptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,&flg);CHKERRQ(ierr); 1321 1322 if (mbs*bs!=m || nbs*bs!=n) { 1323 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Number rows, cols must be divisible by blocksize"); 1324 } 1325 1326 *A = 0; 1327 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,"Mat",comm,MatDestroy,MatView); 1328 PLogObjectCreate(B); 1329 B->data = (void *) (b = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(b); 1330 PetscMemzero(b,sizeof(Mat_SeqBAIJ)); 1331 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 1332 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg); CHKERRQ(ierr); 1333 if (!flg) { 1334 switch (bs) { 1335 case 1: 1336 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 1337 B->ops->solve = MatSolve_SeqBAIJ_1; 1338 B->ops->mult = MatMult_SeqBAIJ_1; 1339 B->ops->multadd = MatMultAdd_SeqBAIJ_1; 1340 break; 1341 case 2: 1342 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 1343 B->ops->solve = MatSolve_SeqBAIJ_2; 1344 B->ops->mult = MatMult_SeqBAIJ_2; 1345 B->ops->multadd = MatMultAdd_SeqBAIJ_2; 1346 break; 1347 case 3: 1348 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 1349 B->ops->solve = MatSolve_SeqBAIJ_3; 1350 B->ops->mult = MatMult_SeqBAIJ_3; 1351 B->ops->multadd = MatMultAdd_SeqBAIJ_3; 1352 break; 1353 case 4: 1354 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 1355 B->ops->solve = MatSolve_SeqBAIJ_4; 1356 B->ops->mult = MatMult_SeqBAIJ_4; 1357 B->ops->multadd = MatMultAdd_SeqBAIJ_4; 1358 break; 1359 case 5: 1360 B->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 1361 B->ops->solve = MatSolve_SeqBAIJ_5; 1362 B->ops->mult = MatMult_SeqBAIJ_5; 1363 B->ops->multadd = MatMultAdd_SeqBAIJ_5; 1364 break; 1365 case 7: 1366 B->ops->mult = MatMult_SeqBAIJ_7; 1367 B->ops->solve = MatSolve_SeqBAIJ_7; 1368 B->ops->multadd = MatMultAdd_SeqBAIJ_7; 1369 break; 1370 } 1371 } 1372 B->ops->destroy = MatDestroy_SeqBAIJ; 1373 B->ops->view = MatView_SeqBAIJ; 1374 B->factor = 0; 1375 B->lupivotthreshold = 1.0; 1376 B->mapping = 0; 1377 b->row = 0; 1378 b->col = 0; 1379 b->icol = 0; 1380 b->reallocs = 0; 1381 1382 b->m = m; B->m = m; B->M = m; 1383 b->n = n; B->n = n; B->N = n; 1384 1385 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1386 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1387 1388 b->mbs = mbs; 1389 b->nbs = nbs; 1390 b->imax = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(b->imax); 1391 if (nnz == PETSC_NULL) { 1392 if (nz == PETSC_DEFAULT) nz = 5; 1393 else if (nz <= 0) nz = 1; 1394 for ( i=0; i<mbs; i++ ) b->imax[i] = nz; 1395 nz = nz*mbs; 1396 } else { 1397 nz = 0; 1398 for ( i=0; i<mbs; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1399 } 1400 1401 /* allocate the matrix space */ 1402 len = nz*sizeof(int) + nz*bs2*sizeof(MatScalar) + (b->m+1)*sizeof(int); 1403 b->a = (MatScalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1404 PetscMemzero(b->a,nz*bs2*sizeof(MatScalar)); 1405 b->j = (int *) (b->a + nz*bs2); 1406 PetscMemzero(b->j,nz*sizeof(int)); 1407 b->i = b->j + nz; 1408 b->singlemalloc = 1; 1409 1410 b->i[0] = 0; 1411 for (i=1; i<mbs+1; i++) { 1412 b->i[i] = b->i[i-1] + b->imax[i-1]; 1413 } 1414 1415 /* b->ilen will count nonzeros in each block row so far. */ 1416 b->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); 1417 PLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1418 for ( i=0; i<mbs; i++ ) { b->ilen[i] = 0;} 1419 1420 b->bs = bs; 1421 b->bs2 = bs2; 1422 b->mbs = mbs; 1423 b->nz = 0; 1424 b->maxnz = nz*bs2; 1425 b->sorted = 0; 1426 b->roworiented = 1; 1427 b->nonew = 0; 1428 b->diag = 0; 1429 b->solve_work = 0; 1430 b->mult_work = 0; 1431 b->spptr = 0; 1432 B->info.nz_unneeded = (double)b->maxnz; 1433 1434 *A = B; 1435 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1436 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1437 1438 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqBAIJSetColumnIndices_C", 1439 "MatSeqBAIJSetColumnIndices_SeqBAIJ", 1440 (void*)MatSeqBAIJSetColumnIndices_SeqBAIJ);CHKERRQ(ierr); 1441 PetscFunctionReturn(0); 1442 } 1443 1444 #undef __FUNC__ 1445 #define __FUNC__ "MatDuplicate_SeqBAIJ" 1446 int MatDuplicate_SeqBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1447 { 1448 Mat C; 1449 Mat_SeqBAIJ *c,*a = (Mat_SeqBAIJ *) A->data; 1450 int i,len, mbs = a->mbs,nz = a->nz,bs2 =a->bs2,ierr; 1451 1452 PetscFunctionBegin; 1453 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,0,"Corrupt matrix"); 1454 1455 *B = 0; 1456 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQBAIJ,"Mat",A->comm,MatDestroy,MatView); 1457 PLogObjectCreate(C); 1458 C->data = (void *) (c = PetscNew(Mat_SeqBAIJ)); CHKPTRQ(c); 1459 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 1460 C->ops->destroy = MatDestroy_SeqBAIJ; 1461 C->ops->view = MatView_SeqBAIJ; 1462 C->factor = A->factor; 1463 c->row = 0; 1464 c->col = 0; 1465 C->assembled = PETSC_TRUE; 1466 1467 c->m = C->m = a->m; 1468 c->n = C->n = a->n; 1469 C->M = a->m; 1470 C->N = a->n; 1471 1472 c->bs = a->bs; 1473 c->bs2 = a->bs2; 1474 c->mbs = a->mbs; 1475 c->nbs = a->nbs; 1476 1477 c->imax = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->imax); 1478 c->ilen = (int *) PetscMalloc((mbs+1)*sizeof(int)); CHKPTRQ(c->ilen); 1479 for ( i=0; i<mbs; i++ ) { 1480 c->imax[i] = a->imax[i]; 1481 c->ilen[i] = a->ilen[i]; 1482 } 1483 1484 /* allocate the matrix space */ 1485 c->singlemalloc = 1; 1486 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1487 c->a = (MatScalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1488 c->j = (int *) (c->a + nz*bs2); 1489 c->i = c->j + nz; 1490 PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int)); 1491 if (mbs > 0) { 1492 PetscMemcpy(c->j,a->j,nz*sizeof(int)); 1493 if (cpvalues == MAT_COPY_VALUES) { 1494 PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar)); 1495 } else { 1496 PetscMemzero(c->a,bs2*nz*sizeof(MatScalar)); 1497 } 1498 } 1499 1500 PLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqBAIJ)); 1501 c->sorted = a->sorted; 1502 c->roworiented = a->roworiented; 1503 c->nonew = a->nonew; 1504 1505 if (a->diag) { 1506 c->diag = (int *) PetscMalloc( (mbs+1)*sizeof(int) ); CHKPTRQ(c->diag); 1507 PLogObjectMemory(C,(mbs+1)*sizeof(int)); 1508 for ( i=0; i<mbs; i++ ) { 1509 c->diag[i] = a->diag[i]; 1510 } 1511 } 1512 else c->diag = 0; 1513 c->nz = a->nz; 1514 c->maxnz = a->maxnz; 1515 c->solve_work = 0; 1516 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1517 c->mult_work = 0; 1518 *B = C; 1519 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1520 PetscFunctionReturn(0); 1521 } 1522 1523 #undef __FUNC__ 1524 #define __FUNC__ "MatLoad_SeqBAIJ" 1525 int MatLoad_SeqBAIJ(Viewer viewer,MatType type,Mat *A) 1526 { 1527 Mat_SeqBAIJ *a; 1528 Mat B; 1529 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1,flg; 1530 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,maskcount; 1531 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1532 int *masked, nmask,tmp,bs2,ishift; 1533 Scalar *aa; 1534 MPI_Comm comm = ((PetscObject) viewer)->comm; 1535 1536 PetscFunctionBegin; 1537 ierr = OptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,&flg);CHKERRQ(ierr); 1538 bs2 = bs*bs; 1539 1540 MPI_Comm_size(comm,&size); 1541 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 1542 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1543 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 1544 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not Mat object"); 1545 M = header[1]; N = header[2]; nz = header[3]; 1546 1547 if (header[3] < 0) { 1548 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format, cannot load as SeqBAIJ"); 1549 } 1550 1551 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Can only do square matrices"); 1552 1553 /* 1554 This code adds extra rows to make sure the number of rows is 1555 divisible by the blocksize 1556 */ 1557 mbs = M/bs; 1558 extra_rows = bs - M + bs*(mbs); 1559 if (extra_rows == bs) extra_rows = 0; 1560 else mbs++; 1561 if (extra_rows) { 1562 PLogInfo(0,"MatLoad_SeqBAIJ:Padding loaded matrix to match blocksize\n"); 1563 } 1564 1565 /* read in row lengths */ 1566 rowlengths = (int*) PetscMalloc((M+extra_rows)*sizeof(int));CHKPTRQ(rowlengths); 1567 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 1568 for ( i=0; i<extra_rows; i++ ) rowlengths[M+i] = 1; 1569 1570 /* read in column indices */ 1571 jj = (int*) PetscMalloc( (nz+extra_rows)*sizeof(int) ); CHKPTRQ(jj); 1572 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT); CHKERRQ(ierr); 1573 for ( i=0; i<extra_rows; i++ ) jj[nz+i] = M+i; 1574 1575 /* loop over row lengths determining block row lengths */ 1576 browlengths = (int *) PetscMalloc(mbs*sizeof(int));CHKPTRQ(browlengths); 1577 PetscMemzero(browlengths,mbs*sizeof(int)); 1578 mask = (int *) PetscMalloc( 2*mbs*sizeof(int) ); CHKPTRQ(mask); 1579 PetscMemzero(mask,mbs*sizeof(int)); 1580 masked = mask + mbs; 1581 rowcount = 0; nzcount = 0; 1582 for ( i=0; i<mbs; i++ ) { 1583 nmask = 0; 1584 for ( j=0; j<bs; j++ ) { 1585 kmax = rowlengths[rowcount]; 1586 for ( k=0; k<kmax; k++ ) { 1587 tmp = jj[nzcount++]/bs; 1588 if (!mask[tmp]) {masked[nmask++] = tmp; mask[tmp] = 1;} 1589 } 1590 rowcount++; 1591 } 1592 browlengths[i] += nmask; 1593 /* zero out the mask elements we set */ 1594 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1595 } 1596 1597 /* create our matrix */ 1598 ierr = MatCreateSeqBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1599 B = *A; 1600 a = (Mat_SeqBAIJ *) B->data; 1601 1602 /* set matrix "i" values */ 1603 a->i[0] = 0; 1604 for ( i=1; i<= mbs; i++ ) { 1605 a->i[i] = a->i[i-1] + browlengths[i-1]; 1606 a->ilen[i-1] = browlengths[i-1]; 1607 } 1608 a->nz = 0; 1609 for ( i=0; i<mbs; i++ ) a->nz += browlengths[i]; 1610 1611 /* read in nonzero values */ 1612 aa = (Scalar *) PetscMalloc((nz+extra_rows)*sizeof(Scalar));CHKPTRQ(aa); 1613 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR); CHKERRQ(ierr); 1614 for ( i=0; i<extra_rows; i++ ) aa[nz+i] = 1.0; 1615 1616 /* set "a" and "j" values into matrix */ 1617 nzcount = 0; jcount = 0; 1618 for ( i=0; i<mbs; i++ ) { 1619 nzcountb = nzcount; 1620 nmask = 0; 1621 for ( j=0; j<bs; j++ ) { 1622 kmax = rowlengths[i*bs+j]; 1623 for ( k=0; k<kmax; k++ ) { 1624 tmp = jj[nzcount++]/bs; 1625 if (!mask[tmp]) { masked[nmask++] = tmp; mask[tmp] = 1;} 1626 } 1627 rowcount++; 1628 } 1629 /* sort the masked values */ 1630 PetscSortInt(nmask,masked); 1631 1632 /* set "j" values into matrix */ 1633 maskcount = 1; 1634 for ( j=0; j<nmask; j++ ) { 1635 a->j[jcount++] = masked[j]; 1636 mask[masked[j]] = maskcount++; 1637 } 1638 /* set "a" values into matrix */ 1639 ishift = bs2*a->i[i]; 1640 for ( j=0; j<bs; j++ ) { 1641 kmax = rowlengths[i*bs+j]; 1642 for ( k=0; k<kmax; k++ ) { 1643 tmp = jj[nzcountb]/bs ; 1644 block = mask[tmp] - 1; 1645 point = jj[nzcountb] - bs*tmp; 1646 idx = ishift + bs2*block + j + bs*point; 1647 a->a[idx] = aa[nzcountb++]; 1648 } 1649 } 1650 /* zero out the mask elements we set */ 1651 for ( j=0; j<nmask; j++ ) mask[masked[j]] = 0; 1652 } 1653 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Bad binary matrix"); 1654 1655 PetscFree(rowlengths); 1656 PetscFree(browlengths); 1657 PetscFree(aa); 1658 PetscFree(jj); 1659 PetscFree(mask); 1660 1661 B->assembled = PETSC_TRUE; 1662 1663 ierr = MatView_Private(B); CHKERRQ(ierr); 1664 PetscFunctionReturn(0); 1665 } 1666 1667 1668 1669