1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: baij.c,v 1.130 1998/04/03 23:15:34 bsmith Exp bsmith $"; 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(Mat A) 106 { 107 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 108 int ierr; 109 110 #if defined(USE_PETSC_LOG) 111 PLogObjectState((PetscObject) A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 112 #endif 113 PetscFree(a->a); 114 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 115 if (a->diag) PetscFree(a->diag); 116 if (a->ilen) PetscFree(a->ilen); 117 if (a->imax) PetscFree(a->imax); 118 if (a->solve_work) PetscFree(a->solve_work); 119 if (a->mult_work) PetscFree(a->mult_work); 120 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 121 PetscFree(a); 122 PLogObjectDestroy(A); 123 PetscHeaderDestroy(A); 124 PetscFunctionReturn(0); 125 } 126 127 #undef __FUNC__ 128 #define __FUNC__ "MatSetOption_SeqBAIJ" 129 int MatSetOption_SeqBAIJ(Mat A,MatOption op) 130 { 131 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 132 133 PetscFunctionBegin; 134 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 135 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 136 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 137 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 138 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 139 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 140 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 141 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 142 else if (op == MAT_ROWS_SORTED || 143 op == MAT_ROWS_UNSORTED || 144 op == MAT_SYMMETRIC || 145 op == MAT_STRUCTURALLY_SYMMETRIC || 146 op == MAT_YES_NEW_DIAGONALS || 147 op == MAT_IGNORE_OFF_PROC_ENTRIES || 148 op == MAT_USE_HASH_TABLE) { 149 PLogInfo(A,"MatSetOption_SeqBAIJ:Option ignored\n"); 150 } else if (op == MAT_NO_NEW_DIAGONALS) { 151 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 152 } else { 153 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 154 } 155 PetscFunctionReturn(0); 156 } 157 158 159 #undef __FUNC__ 160 #define __FUNC__ "MatGetSize_SeqBAIJ" 161 int MatGetSize_SeqBAIJ(Mat A,int *m,int *n) 162 { 163 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 164 165 PetscFunctionBegin; 166 if (m) *m = a->m; 167 if (n) *n = a->n; 168 PetscFunctionReturn(0); 169 } 170 171 #undef __FUNC__ 172 #define __FUNC__ "MatGetOwnershipRange_SeqBAIJ" 173 int MatGetOwnershipRange_SeqBAIJ(Mat A,int *m,int *n) 174 { 175 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 176 177 PetscFunctionBegin; 178 *m = 0; *n = a->m; 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNC__ 183 #define __FUNC__ "MatGetRow_SeqBAIJ" 184 int MatGetRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 185 { 186 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 187 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*idx_i,bs2; 188 Scalar *aa,*v_i,*aa_i; 189 190 PetscFunctionBegin; 191 bs = a->bs; 192 ai = a->i; 193 aj = a->j; 194 aa = a->a; 195 bs2 = a->bs2; 196 197 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 198 199 bn = row/bs; /* Block number */ 200 bp = row % bs; /* Block Position */ 201 M = ai[bn+1] - ai[bn]; 202 *nz = bs*M; 203 204 if (v) { 205 *v = 0; 206 if (*nz) { 207 *v = (Scalar *) PetscMalloc( (*nz)*sizeof(Scalar) ); CHKPTRQ(*v); 208 for ( i=0; i<M; i++ ) { /* for each block in the block row */ 209 v_i = *v + i*bs; 210 aa_i = aa + bs2*(ai[bn] + i); 211 for ( j=bp,k=0; j<bs2; j+=bs,k++ ) {v_i[k] = aa_i[j];} 212 } 213 } 214 } 215 216 if (idx) { 217 *idx = 0; 218 if (*nz) { 219 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 220 for ( i=0; i<M; i++ ) { /* for each block in the block row */ 221 idx_i = *idx + i*bs; 222 itmp = bs*aj[ai[bn] + i]; 223 for ( j=0; j<bs; j++ ) {idx_i[j] = itmp++;} 224 } 225 } 226 } 227 PetscFunctionReturn(0); 228 } 229 230 #undef __FUNC__ 231 #define __FUNC__ "MatRestoreRow_SeqBAIJ" 232 int MatRestoreRow_SeqBAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 233 { 234 PetscFunctionBegin; 235 if (idx) {if (*idx) PetscFree(*idx);} 236 if (v) {if (*v) PetscFree(*v);} 237 PetscFunctionReturn(0); 238 } 239 240 #undef __FUNC__ 241 #define __FUNC__ "MatTranspose_SeqBAIJ" 242 int MatTranspose_SeqBAIJ(Mat A,Mat *B) 243 { 244 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 245 Mat C; 246 int i,j,k,ierr,*aj=a->j,*ai=a->i,bs=a->bs,mbs=a->mbs,nbs=a->nbs,len,*col; 247 int *rows,*cols,bs2=a->bs2; 248 Scalar *array=a->a; 249 250 PetscFunctionBegin; 251 if (B==PETSC_NULL && mbs!=nbs) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Square matrix only for in-place"); 252 col = (int *) PetscMalloc((1+nbs)*sizeof(int)); CHKPTRQ(col); 253 PetscMemzero(col,(1+nbs)*sizeof(int)); 254 255 for ( i=0; i<ai[mbs]; i++ ) col[aj[i]] += 1; 256 ierr = MatCreateSeqBAIJ(A->comm,bs,a->n,a->m,PETSC_NULL,col,&C); CHKERRQ(ierr); 257 PetscFree(col); 258 rows = (int *) PetscMalloc(2*bs*sizeof(int)); CHKPTRQ(rows); 259 cols = rows + bs; 260 for ( i=0; i<mbs; i++ ) { 261 cols[0] = i*bs; 262 for (k=1; k<bs; k++ ) cols[k] = cols[k-1] + 1; 263 len = ai[i+1] - ai[i]; 264 for ( j=0; j<len; j++ ) { 265 rows[0] = (*aj++)*bs; 266 for (k=1; k<bs; k++ ) rows[k] = rows[k-1] + 1; 267 ierr = MatSetValues(C,bs,rows,bs,cols,array,INSERT_VALUES); CHKERRQ(ierr); 268 array += bs2; 269 } 270 } 271 PetscFree(rows); 272 273 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 274 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 275 276 if (B != PETSC_NULL) { 277 *B = C; 278 } else { 279 PetscOps *Abops; 280 struct _MatOps *Aops; 281 282 /* This isn't really an in-place transpose */ 283 PetscFree(a->a); 284 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 285 if (a->diag) PetscFree(a->diag); 286 if (a->ilen) PetscFree(a->ilen); 287 if (a->imax) PetscFree(a->imax); 288 if (a->solve_work) PetscFree(a->solve_work); 289 PetscFree(a); 290 291 /* 292 This is horrible, horrible code. We need to keep the 293 A pointers for the bops and ops but copy everything 294 else from C. 295 */ 296 Abops = A->bops; 297 Aops = A->ops; 298 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 299 A->bops = Abops; 300 A->ops = Aops; 301 302 PetscHeaderDestroy(C); 303 } 304 PetscFunctionReturn(0); 305 } 306 307 308 309 310 #undef __FUNC__ 311 #define __FUNC__ "MatView_SeqBAIJ_Binary" 312 static int MatView_SeqBAIJ_Binary(Mat A,Viewer viewer) 313 { 314 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 315 int i, fd, *col_lens, ierr, bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 316 Scalar *aa; 317 318 PetscFunctionBegin; 319 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 320 col_lens = (int *) PetscMalloc((4+a->m)*sizeof(int));CHKPTRQ(col_lens); 321 col_lens[0] = MAT_COOKIE; 322 323 col_lens[1] = a->m; 324 col_lens[2] = a->n; 325 col_lens[3] = a->nz*bs2; 326 327 /* store lengths of each row and write (including header) to file */ 328 count = 0; 329 for ( i=0; i<a->mbs; i++ ) { 330 for ( j=0; j<bs; j++ ) { 331 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 332 } 333 } 334 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1); CHKERRQ(ierr); 335 PetscFree(col_lens); 336 337 /* store column indices (zero start index) */ 338 jj = (int *) PetscMalloc( (a->nz+1)*bs2*sizeof(int) ); CHKPTRQ(jj); 339 count = 0; 340 for ( i=0; i<a->mbs; i++ ) { 341 for ( j=0; j<bs; j++ ) { 342 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 343 for ( l=0; l<bs; l++ ) { 344 jj[count++] = bs*a->j[k] + l; 345 } 346 } 347 } 348 } 349 ierr = PetscBinaryWrite(fd,jj,bs2*a->nz,PETSC_INT,0); CHKERRQ(ierr); 350 PetscFree(jj); 351 352 /* store nonzero values */ 353 aa = (Scalar *) PetscMalloc((a->nz+1)*bs2*sizeof(Scalar)); CHKPTRQ(aa); 354 count = 0; 355 for ( i=0; i<a->mbs; i++ ) { 356 for ( j=0; j<bs; j++ ) { 357 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 358 for ( l=0; l<bs; l++ ) { 359 aa[count++] = a->a[bs2*k + l*bs + j]; 360 } 361 } 362 } 363 } 364 ierr = PetscBinaryWrite(fd,aa,bs2*a->nz,PETSC_SCALAR,0); CHKERRQ(ierr); 365 PetscFree(aa); 366 PetscFunctionReturn(0); 367 } 368 369 #undef __FUNC__ 370 #define __FUNC__ "MatView_SeqBAIJ_ASCII" 371 static int MatView_SeqBAIJ_ASCII(Mat A,Viewer viewer) 372 { 373 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data; 374 int ierr, i,j,format,bs = a->bs,k,l,bs2=a->bs2; 375 FILE *fd; 376 char *outputname; 377 378 PetscFunctionBegin; 379 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 380 ierr = ViewerFileGetOutputname_Private(viewer,&outputname);CHKERRQ(ierr); 381 ierr = ViewerGetFormat(viewer,&format); 382 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 383 fprintf(fd," block size is %d\n",bs); 384 } 385 else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 386 SETERRQ(PETSC_ERR_SUP,0,"Matlab format not supported"); 387 } 388 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 389 for ( i=0; i<a->mbs; i++ ) { 390 for ( j=0; j<bs; j++ ) { 391 fprintf(fd,"row %d:",i*bs+j); 392 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 393 for ( l=0; l<bs; l++ ) { 394 #if defined(USE_PETSC_COMPLEX) 395 if (imag(a->a[bs2*k + l*bs + j]) > 0.0 && real(a->a[bs2*k + l*bs + j]) != 0.0) 396 fprintf(fd," %d %g + %g i",bs*a->j[k]+l, 397 real(a->a[bs2*k + l*bs + j]),imag(a->a[bs2*k + l*bs + j])); 398 else if (imag(a->a[bs2*k + l*bs + j]) < 0.0 && real(a->a[bs2*k + l*bs + j]) != 0.0) 399 fprintf(fd," %d %g - %g i",bs*a->j[k]+l, 400 real(a->a[bs2*k + l*bs + j]),-imag(a->a[bs2*k + l*bs + j])); 401 else if (real(a->a[bs2*k + l*bs + j]) != 0.0) 402 fprintf(fd," %d %g ",bs*a->j[k]+l,real(a->a[bs2*k + l*bs + j])); 403 #else 404 if (a->a[bs2*k + l*bs + j] != 0.0) 405 fprintf(fd," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]); 406 #endif 407 } 408 } 409 fprintf(fd,"\n"); 410 } 411 } 412 } 413 else { 414 for ( i=0; i<a->mbs; i++ ) { 415 for ( j=0; j<bs; j++ ) { 416 fprintf(fd,"row %d:",i*bs+j); 417 for ( k=a->i[i]; k<a->i[i+1]; k++ ) { 418 for ( l=0; l<bs; l++ ) { 419 #if defined(USE_PETSC_COMPLEX) 420 if (imag(a->a[bs2*k + l*bs + j]) > 0.0) { 421 fprintf(fd," %d %g + %g i",bs*a->j[k]+l, 422 real(a->a[bs2*k + l*bs + j]),imag(a->a[bs2*k + l*bs + j])); 423 } 424 else if (imag(a->a[bs2*k + l*bs + j]) < 0.0) { 425 fprintf(fd," %d %g - %g i",bs*a->j[k]+l, 426 real(a->a[bs2*k + l*bs + j]),-imag(a->a[bs2*k + l*bs + j])); 427 } 428 else { 429 fprintf(fd," %d %g ",bs*a->j[k]+l,real(a->a[bs2*k + l*bs + j])); 430 } 431 #else 432 fprintf(fd," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]); 433 #endif 434 } 435 } 436 fprintf(fd,"\n"); 437 } 438 } 439 } 440 fflush(fd); 441 PetscFunctionReturn(0); 442 } 443 444 #undef __FUNC__ 445 #define __FUNC__ "MatView_SeqBAIJ_Draw" 446 static int MatView_SeqBAIJ_Draw(Mat A,Viewer viewer) 447 { 448 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *) A->data; 449 int row,ierr,i,j,k,l,mbs=a->mbs,pause,color,bs=a->bs,bs2=a->bs2; 450 double xl,yl,xr,yr,w,h,xc,yc,scale = 1.0,x_l,x_r,y_l,y_r; 451 Scalar *aa; 452 Draw draw; 453 DrawButton button; 454 PetscTruth isnull; 455 456 PetscFunctionBegin; 457 ierr = ViewerDrawGetDraw(viewer,&draw);CHKERRQ(ierr); 458 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 459 460 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 461 xr += w; yr += h; xl = -w; yl = -h; 462 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 463 /* loop over matrix elements drawing boxes */ 464 color = DRAW_BLUE; 465 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 466 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 467 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 468 x_l = a->j[j]*bs; x_r = x_l + 1.0; 469 aa = a->a + j*bs2; 470 for ( k=0; k<bs; k++ ) { 471 for ( l=0; l<bs; l++ ) { 472 if (PetscReal(*aa++) >= 0.) continue; 473 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 474 } 475 } 476 } 477 } 478 color = DRAW_CYAN; 479 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 480 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 481 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 482 x_l = a->j[j]*bs; x_r = x_l + 1.0; 483 aa = a->a + j*bs2; 484 for ( k=0; k<bs; k++ ) { 485 for ( l=0; l<bs; l++ ) { 486 if (PetscReal(*aa++) != 0.) continue; 487 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 488 } 489 } 490 } 491 } 492 493 color = DRAW_RED; 494 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 495 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 496 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 497 x_l = a->j[j]*bs; x_r = x_l + 1.0; 498 aa = a->a + j*bs2; 499 for ( k=0; k<bs; k++ ) { 500 for ( l=0; l<bs; l++ ) { 501 if (PetscReal(*aa++) <= 0.) continue; 502 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 503 } 504 } 505 } 506 } 507 508 DrawSynchronizedFlush(draw); 509 DrawGetPause(draw,&pause); 510 if (pause >= 0) { PetscSleep(pause); PetscFunctionReturn(0);} 511 512 /* allow the matrix to zoom or shrink */ 513 ierr = DrawSynchronizedGetMouseButton(draw,&button,&xc,&yc,0,0); 514 while (button != BUTTON_RIGHT) { 515 DrawSynchronizedClear(draw); 516 if (button == BUTTON_LEFT) scale = .5; 517 else if (button == BUTTON_CENTER) scale = 2.; 518 xl = scale*(xl + w - xc) + xc - w*scale; 519 xr = scale*(xr - w - xc) + xc + w*scale; 520 yl = scale*(yl + h - yc) + yc - h*scale; 521 yr = scale*(yr - h - yc) + yc + h*scale; 522 w *= scale; h *= scale; 523 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 524 color = DRAW_BLUE; 525 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 526 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 527 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 528 x_l = a->j[j]*bs; x_r = x_l + 1.0; 529 aa = a->a + j*bs2; 530 for ( k=0; k<bs; k++ ) { 531 for ( l=0; l<bs; l++ ) { 532 if (PetscReal(*aa++) >= 0.) continue; 533 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 534 } 535 } 536 } 537 } 538 color = DRAW_CYAN; 539 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 540 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 541 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 542 x_l = a->j[j]*bs; x_r = x_l + 1.0; 543 aa = a->a + j*bs2; 544 for ( k=0; k<bs; k++ ) { 545 for ( l=0; l<bs; l++ ) { 546 if (PetscReal(*aa++) != 0.) continue; 547 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 548 } 549 } 550 } 551 } 552 553 color = DRAW_RED; 554 for ( i=0,row=0; i<mbs; i++,row+=bs ) { 555 for ( j=a->i[i]; j<a->i[i+1]; j++ ) { 556 y_l = a->m - row - 1.0; y_r = y_l + 1.0; 557 x_l = a->j[j]*bs; x_r = x_l + 1.0; 558 aa = a->a + j*bs2; 559 for ( k=0; k<bs; k++ ) { 560 for ( l=0; l<bs; l++ ) { 561 if (PetscReal(*aa++) <= 0.) continue; 562 DrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color); 563 } 564 } 565 } 566 } 567 ierr = DrawSynchronizedGetMouseButton(draw,&button,&xc,&yc,0,0); 568 } 569 PetscFunctionReturn(0); 570 } 571 572 #undef __FUNC__ 573 #define __FUNC__ "MatView_SeqBAIJ" 574 int MatView_SeqBAIJ(Mat A,Viewer viewer) 575 { 576 ViewerType vtype; 577 int ierr; 578 579 PetscFunctionBegin; 580 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 581 if (vtype == MATLAB_VIEWER) { 582 SETERRQ(PETSC_ERR_SUP,0,"Matlab viewer not supported"); 583 } else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER){ 584 ierr = MatView_SeqBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 585 } else if (vtype == BINARY_FILE_VIEWER) { 586 ierr = MatView_SeqBAIJ_Binary(A,viewer);CHKERRQ(ierr); 587 } else if (vtype == DRAW_VIEWER) { 588 ierr = MatView_SeqBAIJ_Draw(A,viewer);CHKERRQ(ierr); 589 } else { 590 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 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->ops->destroy = MatDestroy_SeqBAIJ; 1260 B->ops->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->ops->destroy = MatDestroy_SeqBAIJ; 1340 C->ops->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