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