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