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