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