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