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