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