1 /*$Id: sbaij.c,v 1.62 2001/08/07 03:03:01 balay Exp $*/ 2 3 /* 4 Defines the basic matrix operations for the BAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "src/mat/impls/baij/seq/baij.h" /*I "petscmat.h" I*/ 8 #include "src/vec/vecimpl.h" 9 #include "src/inline/spops.h" 10 #include "src/mat/impls/sbaij/seq/sbaij.h" 11 12 #define CHUNKSIZE 10 13 14 /* 15 Checks for missing diagonals 16 */ 17 #undef __FUNCT__ 18 #define __FUNCT__ "MatMissingDiagonal_SeqSBAIJ" 19 int MatMissingDiagonal_SeqSBAIJ(Mat A) 20 { 21 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 22 int *diag,*jj = a->j,i,ierr; 23 24 PetscFunctionBegin; 25 ierr = MatMarkDiagonal_SeqSBAIJ(A);CHKERRQ(ierr); 26 diag = a->diag; 27 for (i=0; i<a->mbs; i++) { 28 if (jj[diag[i]] != i) { 29 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 30 } 31 } 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "MatMarkDiagonal_SeqSBAIJ" 37 int MatMarkDiagonal_SeqSBAIJ(Mat A) 38 { 39 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 40 int i,j,*diag,m = a->mbs,ierr; 41 42 PetscFunctionBegin; 43 if (a->diag) PetscFunctionReturn(0); 44 45 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 46 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 47 for (i=0; i<m; i++) { 48 diag[i] = a->i[i+1]; 49 for (j=a->i[i]; j<a->i[i+1]; j++) { 50 if (a->j[j] == i) { 51 diag[i] = j; 52 break; 53 } 54 } 55 } 56 a->diag = diag; 57 PetscFunctionReturn(0); 58 } 59 60 61 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 62 63 #undef __FUNCT__ 64 #define __FUNCT__ "MatGetRowIJ_SeqSBAIJ" 65 static int MatGetRowIJ_SeqSBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 66 { 67 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 68 69 PetscFunctionBegin; 70 if (ia) { 71 SETERRQ(1,"Function not yet written for SBAIJ format, only supports natural ordering"); 72 } 73 *nn = a->mbs; 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "MatRestoreRowIJ_SeqSBAIJ" 79 static int MatRestoreRowIJ_SeqSBAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja,PetscTruth *done) 80 { 81 PetscFunctionBegin; 82 if (!ia) PetscFunctionReturn(0); 83 SETERRQ(1,"Function not yet written for SBAIJ format, only supports natural ordering"); 84 /* PetscFunctionReturn(0); */ 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "MatGetBlockSize_SeqSBAIJ" 89 int MatGetBlockSize_SeqSBAIJ(Mat mat,int *bs) 90 { 91 Mat_SeqSBAIJ *sbaij = (Mat_SeqSBAIJ*)mat->data; 92 93 PetscFunctionBegin; 94 *bs = sbaij->bs; 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "MatDestroy_SeqSBAIJ" 100 int MatDestroy_SeqSBAIJ(Mat A) 101 { 102 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 103 int ierr; 104 105 PetscFunctionBegin; 106 #if defined(PETSC_USE_LOG) 107 PetscLogObjectState((PetscObject)A,"Rows=%d, s_NZ=%d",A->m,a->s_nz); 108 #endif 109 ierr = PetscFree(a->a);CHKERRQ(ierr); 110 if (!a->singlemalloc) { 111 ierr = PetscFree(a->i);CHKERRQ(ierr); 112 ierr = PetscFree(a->j);CHKERRQ(ierr); 113 } 114 if (a->row) { 115 ierr = ISDestroy(a->row);CHKERRQ(ierr); 116 } 117 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 118 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 119 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 120 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 121 if (a->mult_work) {ierr = PetscFree(a->mult_work);CHKERRQ(ierr);} 122 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 123 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 124 125 if (a->inew){ 126 ierr = PetscFree(a->inew);CHKERRQ(ierr); 127 a->inew = 0; 128 } 129 ierr = PetscFree(a);CHKERRQ(ierr); 130 PetscFunctionReturn(0); 131 } 132 133 #undef __FUNCT__ 134 #define __FUNCT__ "MatSetOption_SeqSBAIJ" 135 int MatSetOption_SeqSBAIJ(Mat A,MatOption op) 136 { 137 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 138 139 PetscFunctionBegin; 140 switch (op) { 141 case MAT_ROW_ORIENTED: 142 a->roworiented = PETSC_TRUE; 143 break; 144 case MAT_COLUMN_ORIENTED: 145 a->roworiented = PETSC_FALSE; 146 break; 147 case MAT_COLUMNS_SORTED: 148 a->sorted = PETSC_TRUE; 149 break; 150 case MAT_COLUMNS_UNSORTED: 151 a->sorted = PETSC_FALSE; 152 break; 153 case MAT_KEEP_ZEROED_ROWS: 154 a->keepzeroedrows = PETSC_TRUE; 155 break; 156 case MAT_NO_NEW_NONZERO_LOCATIONS: 157 a->nonew = 1; 158 break; 159 case MAT_NEW_NONZERO_LOCATION_ERR: 160 a->nonew = -1; 161 break; 162 case MAT_NEW_NONZERO_ALLOCATION_ERR: 163 a->nonew = -2; 164 break; 165 case MAT_YES_NEW_NONZERO_LOCATIONS: 166 a->nonew = 0; 167 break; 168 case MAT_ROWS_SORTED: 169 case MAT_ROWS_UNSORTED: 170 case MAT_YES_NEW_DIAGONALS: 171 case MAT_IGNORE_OFF_PROC_ENTRIES: 172 case MAT_USE_HASH_TABLE: 173 case MAT_USE_SINGLE_PRECISION_SOLVES: 174 PetscLogInfo(A,"MatSetOption_SeqSBAIJ:Option ignored\n"); 175 break; 176 case MAT_NO_NEW_DIAGONALS: 177 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 178 break; 179 default: 180 SETERRQ(PETSC_ERR_SUP,"unknown option"); 181 break; 182 } 183 PetscFunctionReturn(0); 184 } 185 186 #undef __FUNCT__ 187 #define __FUNCT__ "MatGetRow_SeqSBAIJ" 188 int MatGetRow_SeqSBAIJ(Mat A,int row,int *ncols,int **cols,PetscScalar **v) 189 { 190 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 191 int itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2,ierr; 192 MatScalar *aa,*aa_i; 193 PetscScalar *v_i; 194 195 PetscFunctionBegin; 196 bs = a->bs; 197 ai = a->i; 198 aj = a->j; 199 aa = a->a; 200 bs2 = a->bs2; 201 202 if (row < 0 || row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row out of range"); 203 204 bn = row/bs; /* Block number */ 205 bp = row % bs; /* Block position */ 206 M = ai[bn+1] - ai[bn]; 207 *ncols = bs*M; 208 209 if (v) { 210 *v = 0; 211 if (*ncols) { 212 ierr = PetscMalloc((*ncols+row)*sizeof(PetscScalar),v);CHKERRQ(ierr); 213 for (i=0; i<M; i++) { /* for each block in the block row */ 214 v_i = *v + i*bs; 215 aa_i = aa + bs2*(ai[bn] + i); 216 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 217 } 218 } 219 } 220 221 if (cols) { 222 *cols = 0; 223 if (*ncols) { 224 ierr = PetscMalloc((*ncols+row)*sizeof(int),cols);CHKERRQ(ierr); 225 for (i=0; i<M; i++) { /* for each block in the block row */ 226 cols_i = *cols + i*bs; 227 itmp = bs*aj[ai[bn] + i]; 228 for (j=0; j<bs; j++) {cols_i[j] = itmp++;} 229 } 230 } 231 } 232 233 /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */ 234 /* this segment is currently removed, so only entries in the upper triangle are obtained */ 235 #ifdef column_search 236 v_i = *v + M*bs; 237 cols_i = *cols + M*bs; 238 for (i=0; i<bn; i++){ /* for each block row */ 239 M = ai[i+1] - ai[i]; 240 for (j=0; j<M; j++){ 241 itmp = aj[ai[i] + j]; /* block column value */ 242 if (itmp == bn){ 243 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 244 for (k=0; k<bs; k++) { 245 *cols_i++ = i*bs+k; 246 *v_i++ = aa_i[k]; 247 } 248 *ncols += bs; 249 break; 250 } 251 } 252 } 253 #endif 254 255 PetscFunctionReturn(0); 256 } 257 258 #undef __FUNCT__ 259 #define __FUNCT__ "MatRestoreRow_SeqSBAIJ" 260 int MatRestoreRow_SeqSBAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 261 { 262 int ierr; 263 264 PetscFunctionBegin; 265 if (idx) {if (*idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 266 if (v) {if (*v) {ierr = PetscFree(*v);CHKERRQ(ierr);}} 267 PetscFunctionReturn(0); 268 } 269 270 #undef __FUNCT__ 271 #define __FUNCT__ "MatTranspose_SeqSBAIJ" 272 int MatTranspose_SeqSBAIJ(Mat A,Mat *B) 273 { 274 PetscFunctionBegin; 275 SETERRQ(1,"Matrix is symmetric. MatTranspose() should not be called"); 276 /* PetscFunctionReturn(0); */ 277 } 278 279 #undef __FUNCT__ 280 #define __FUNCT__ "MatView_SeqSBAIJ_Binary" 281 static int MatView_SeqSBAIJ_Binary(Mat A,PetscViewer viewer) 282 { 283 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 284 int i,fd,*col_lens,ierr,bs = a->bs,count,*jj,j,k,l,bs2=a->bs2; 285 PetscScalar *aa; 286 FILE *file; 287 288 PetscFunctionBegin; 289 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 290 ierr = PetscMalloc((4+A->m)*sizeof(int),&col_lens);CHKERRQ(ierr); 291 col_lens[0] = MAT_FILE_COOKIE; 292 293 col_lens[1] = A->m; 294 col_lens[2] = A->m; 295 col_lens[3] = a->s_nz*bs2; 296 297 /* store lengths of each row and write (including header) to file */ 298 count = 0; 299 for (i=0; i<a->mbs; i++) { 300 for (j=0; j<bs; j++) { 301 col_lens[4+count++] = bs*(a->i[i+1] - a->i[i]); 302 } 303 } 304 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,1);CHKERRQ(ierr); 305 ierr = PetscFree(col_lens);CHKERRQ(ierr); 306 307 /* store column indices (zero start index) */ 308 ierr = PetscMalloc((a->s_nz+1)*bs2*sizeof(int),&jj);CHKERRQ(ierr); 309 count = 0; 310 for (i=0; i<a->mbs; i++) { 311 for (j=0; j<bs; j++) { 312 for (k=a->i[i]; k<a->i[i+1]; k++) { 313 for (l=0; l<bs; l++) { 314 jj[count++] = bs*a->j[k] + l; 315 } 316 } 317 } 318 } 319 ierr = PetscBinaryWrite(fd,jj,bs2*a->s_nz,PETSC_INT,0);CHKERRQ(ierr); 320 ierr = PetscFree(jj);CHKERRQ(ierr); 321 322 /* store nonzero values */ 323 ierr = PetscMalloc((a->s_nz+1)*bs2*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 324 count = 0; 325 for (i=0; i<a->mbs; i++) { 326 for (j=0; j<bs; j++) { 327 for (k=a->i[i]; k<a->i[i+1]; k++) { 328 for (l=0; l<bs; l++) { 329 aa[count++] = a->a[bs2*k + l*bs + j]; 330 } 331 } 332 } 333 } 334 ierr = PetscBinaryWrite(fd,aa,bs2*a->s_nz,PETSC_SCALAR,0);CHKERRQ(ierr); 335 ierr = PetscFree(aa);CHKERRQ(ierr); 336 337 ierr = PetscViewerBinaryGetInfoPointer(viewer,&file);CHKERRQ(ierr); 338 if (file) { 339 fprintf(file,"-matload_block_size %d\n",a->bs); 340 } 341 PetscFunctionReturn(0); 342 } 343 344 #undef __FUNCT__ 345 #define __FUNCT__ "MatView_SeqSBAIJ_ASCII" 346 static int MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 347 { 348 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 349 int ierr,i,j,bs = a->bs,k,l,bs2=a->bs2; 350 char *name; 351 PetscViewerFormat format; 352 353 PetscFunctionBegin; 354 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 355 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 356 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_LONG) { 357 ierr = PetscViewerASCIIPrintf(viewer," block size is %d\n",bs);CHKERRQ(ierr); 358 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 359 SETERRQ(PETSC_ERR_SUP,"Matlab format not supported"); 360 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 361 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 362 for (i=0; i<a->mbs; i++) { 363 for (j=0; j<bs; j++) { 364 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 365 for (k=a->i[i]; k<a->i[i+1]; k++) { 366 for (l=0; l<bs; l++) { 367 #if defined(PETSC_USE_COMPLEX) 368 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 369 ierr = PetscViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 370 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 371 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 372 ierr = PetscViewerASCIIPrintf(viewer," %d %g - %g i",bs*a->j[k]+l, 373 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 374 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 375 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 376 } 377 #else 378 if (a->a[bs2*k + l*bs + j] != 0.0) { 379 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 380 } 381 #endif 382 } 383 } 384 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 385 } 386 } 387 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 388 } else { 389 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 390 for (i=0; i<a->mbs; i++) { 391 for (j=0; j<bs; j++) { 392 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i*bs+j);CHKERRQ(ierr); 393 for (k=a->i[i]; k<a->i[i+1]; k++) { 394 for (l=0; l<bs; l++) { 395 #if defined(PETSC_USE_COMPLEX) 396 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 397 ierr = PetscViewerASCIIPrintf(viewer," %d %g + %g i",bs*a->j[k]+l, 398 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 399 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 400 ierr = PetscViewerASCIIPrintf(viewer," %d %g - %g i",bs*a->j[k]+l, 401 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 402 } else { 403 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 404 } 405 #else 406 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 407 #endif 408 } 409 } 410 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 411 } 412 } 413 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 414 } 415 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 416 PetscFunctionReturn(0); 417 } 418 419 #undef __FUNCT__ 420 #define __FUNCT__ "MatView_SeqSBAIJ_Draw_Zoom" 421 static int MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 422 { 423 Mat A = (Mat) Aa; 424 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 425 int row,ierr,i,j,k,l,mbs=a->mbs,color,bs=a->bs,bs2=a->bs2,rank; 426 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 427 MatScalar *aa; 428 MPI_Comm comm; 429 PetscViewer viewer; 430 431 PetscFunctionBegin; 432 /* 433 This is nasty. If this is called from an originally parallel matrix 434 then all processes call this,but only the first has the matrix so the 435 rest should return immediately. 436 */ 437 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 438 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 439 if (rank) PetscFunctionReturn(0); 440 441 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 442 443 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 444 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 445 446 /* loop over matrix elements drawing boxes */ 447 color = PETSC_DRAW_BLUE; 448 for (i=0,row=0; i<mbs; i++,row+=bs) { 449 for (j=a->i[i]; j<a->i[i+1]; j++) { 450 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 451 x_l = a->j[j]*bs; x_r = x_l + 1.0; 452 aa = a->a + j*bs2; 453 for (k=0; k<bs; k++) { 454 for (l=0; l<bs; l++) { 455 if (PetscRealPart(*aa++) >= 0.) continue; 456 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 457 } 458 } 459 } 460 } 461 color = PETSC_DRAW_CYAN; 462 for (i=0,row=0; i<mbs; i++,row+=bs) { 463 for (j=a->i[i]; j<a->i[i+1]; j++) { 464 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 465 x_l = a->j[j]*bs; x_r = x_l + 1.0; 466 aa = a->a + j*bs2; 467 for (k=0; k<bs; k++) { 468 for (l=0; l<bs; l++) { 469 if (PetscRealPart(*aa++) != 0.) continue; 470 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 471 } 472 } 473 } 474 } 475 476 color = PETSC_DRAW_RED; 477 for (i=0,row=0; i<mbs; i++,row+=bs) { 478 for (j=a->i[i]; j<a->i[i+1]; j++) { 479 y_l = A->m - row - 1.0; y_r = y_l + 1.0; 480 x_l = a->j[j]*bs; x_r = x_l + 1.0; 481 aa = a->a + j*bs2; 482 for (k=0; k<bs; k++) { 483 for (l=0; l<bs; l++) { 484 if (PetscRealPart(*aa++) <= 0.) continue; 485 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 486 } 487 } 488 } 489 } 490 PetscFunctionReturn(0); 491 } 492 493 #undef __FUNCT__ 494 #define __FUNCT__ "MatView_SeqSBAIJ_Draw" 495 static int MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 496 { 497 int ierr; 498 PetscReal xl,yl,xr,yr,w,h; 499 PetscDraw draw; 500 PetscTruth isnull; 501 502 PetscFunctionBegin; 503 504 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 505 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 506 507 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 508 xr = A->m; yr = A->m; h = yr/10.0; w = xr/10.0; 509 xr += w; yr += h; xl = -w; yl = -h; 510 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 511 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 512 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 513 PetscFunctionReturn(0); 514 } 515 516 #undef __FUNCT__ 517 #define __FUNCT__ "MatView_SeqSBAIJ" 518 int MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 519 { 520 int ierr; 521 PetscTruth issocket,isascii,isbinary,isdraw; 522 523 PetscFunctionBegin; 524 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 525 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 526 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 527 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 528 if (issocket) { 529 SETERRQ(PETSC_ERR_SUP,"Socket viewer not supported"); 530 } else if (isascii){ 531 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 532 } else if (isbinary) { 533 ierr = MatView_SeqSBAIJ_Binary(A,viewer);CHKERRQ(ierr); 534 } else if (isdraw) { 535 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 536 } else { 537 SETERRQ1(1,"Viewer type %s not supported by SeqBAIJ matrices",((PetscObject)viewer)->type_name); 538 } 539 PetscFunctionReturn(0); 540 } 541 542 543 #undef __FUNCT__ 544 #define __FUNCT__ "MatGetValues_SeqSBAIJ" 545 int MatGetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v) 546 { 547 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 548 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 549 int *ai = a->i,*ailen = a->ilen; 550 int brow,bcol,ridx,cidx,bs=a->bs,bs2=a->bs2; 551 MatScalar *ap,*aa = a->a,zero = 0.0; 552 553 PetscFunctionBegin; 554 for (k=0; k<m; k++) { /* loop over rows */ 555 row = im[k]; brow = row/bs; 556 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative row"); 557 if (row >= A->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 558 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 559 nrow = ailen[brow]; 560 for (l=0; l<n; l++) { /* loop over columns */ 561 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative column"); 562 if (in[l] >= A->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 563 col = in[l] ; 564 bcol = col/bs; 565 cidx = col%bs; 566 ridx = row%bs; 567 high = nrow; 568 low = 0; /* assume unsorted */ 569 while (high-low > 5) { 570 t = (low+high)/2; 571 if (rp[t] > bcol) high = t; 572 else low = t; 573 } 574 for (i=low; i<high; i++) { 575 if (rp[i] > bcol) break; 576 if (rp[i] == bcol) { 577 *v++ = ap[bs2*i+bs*cidx+ridx]; 578 goto finished; 579 } 580 } 581 *v++ = zero; 582 finished:; 583 } 584 } 585 PetscFunctionReturn(0); 586 } 587 588 589 #undef __FUNCT__ 590 #define __FUNCT__ "MatSetValuesBlocked_SeqSBAIJ" 591 int MatSetValuesBlocked_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v,InsertMode is) 592 { 593 PetscFunctionBegin; 594 SETERRQ(1,"Function not yet written for SBAIJ format"); 595 /* PetscFunctionReturn(0); */ 596 } 597 598 #undef __FUNCT__ 599 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 600 int MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 601 { 602 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 603 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 604 int m = A->m,*ip,N,*ailen = a->ilen; 605 int mbs = a->mbs,bs2 = a->bs2,rmax = 0,ierr; 606 MatScalar *aa = a->a,*ap; 607 608 PetscFunctionBegin; 609 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 610 611 if (m) rmax = ailen[0]; 612 for (i=1; i<mbs; i++) { 613 /* move each row back by the amount of empty slots (fshift) before it*/ 614 fshift += imax[i-1] - ailen[i-1]; 615 rmax = PetscMax(rmax,ailen[i]); 616 if (fshift) { 617 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 618 N = ailen[i]; 619 for (j=0; j<N; j++) { 620 ip[j-fshift] = ip[j]; 621 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 622 } 623 } 624 ai[i] = ai[i-1] + ailen[i-1]; 625 } 626 if (mbs) { 627 fshift += imax[mbs-1] - ailen[mbs-1]; 628 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 629 } 630 /* reset ilen and imax for each row */ 631 for (i=0; i<mbs; i++) { 632 ailen[i] = imax[i] = ai[i+1] - ai[i]; 633 } 634 a->s_nz = ai[mbs]; 635 636 /* diagonals may have moved, so kill the diagonal pointers */ 637 if (fshift && a->diag) { 638 ierr = PetscFree(a->diag);CHKERRQ(ierr); 639 PetscLogObjectMemory(A,-(m+1)*sizeof(int)); 640 a->diag = 0; 641 } 642 PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Matrix size: %d X %d, block size %d; storage space: %d unneeded, %d used\n", 643 m,A->m,a->bs,fshift*bs2,a->s_nz*bs2); 644 PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Number of mallocs during MatSetValues is %d\n", 645 a->reallocs); 646 PetscLogInfo(A,"MatAssemblyEnd_SeqSBAIJ:Most nonzeros blocks in any row is %d\n",rmax); 647 a->reallocs = 0; 648 A->info.nz_unneeded = (PetscReal)fshift*bs2; 649 650 PetscFunctionReturn(0); 651 } 652 653 /* 654 This function returns an array of flags which indicate the locations of contiguous 655 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 656 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 657 Assume: sizes should be long enough to hold all the values. 658 */ 659 #undef __FUNCT__ 660 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 661 static int MatZeroRows_SeqSBAIJ_Check_Blocks(int idx[],int n,int bs,int sizes[], int *bs_max) 662 { 663 int i,j,k,row; 664 PetscTruth flg; 665 666 PetscFunctionBegin; 667 for (i=0,j=0; i<n; j++) { 668 row = idx[i]; 669 if (row%bs!=0) { /* Not the begining of a block */ 670 sizes[j] = 1; 671 i++; 672 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 673 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 674 i++; 675 } else { /* Begining of the block, so check if the complete block exists */ 676 flg = PETSC_TRUE; 677 for (k=1; k<bs; k++) { 678 if (row+k != idx[i+k]) { /* break in the block */ 679 flg = PETSC_FALSE; 680 break; 681 } 682 } 683 if (flg == PETSC_TRUE) { /* No break in the bs */ 684 sizes[j] = bs; 685 i+= bs; 686 } else { 687 sizes[j] = 1; 688 i++; 689 } 690 } 691 } 692 *bs_max = j; 693 PetscFunctionReturn(0); 694 } 695 696 #undef __FUNCT__ 697 #define __FUNCT__ "MatZeroRows_SeqSBAIJ" 698 int MatZeroRows_SeqSBAIJ(Mat A,IS is,PetscScalar *diag) 699 { 700 Mat_SeqSBAIJ *sbaij=(Mat_SeqSBAIJ*)A->data; 701 int ierr,i,j,k,count,is_n,*is_idx,*rows; 702 int bs=sbaij->bs,bs2=sbaij->bs2,*sizes,row,bs_max; 703 PetscScalar zero = 0.0; 704 MatScalar *aa; 705 706 PetscFunctionBegin; 707 /* Make a copy of the IS and sort it */ 708 ierr = ISGetSize(is,&is_n);CHKERRQ(ierr); 709 ierr = ISGetIndices(is,&is_idx);CHKERRQ(ierr); 710 711 /* allocate memory for rows,sizes */ 712 ierr = PetscMalloc((3*is_n+1)*sizeof(int),&rows);CHKERRQ(ierr); 713 sizes = rows + is_n; 714 715 /* initialize copy IS values to rows, and sort them */ 716 for (i=0; i<is_n; i++) { rows[i] = is_idx[i]; } 717 ierr = PetscSortInt(is_n,rows);CHKERRQ(ierr); 718 if (sbaij->keepzeroedrows) { /* do not change nonzero structure */ 719 for (i=0; i<is_n; i++) { sizes[i] = 1; } /* sizes: size of blocks, = 1 or bs */ 720 bs_max = is_n; /* bs_max: num. of contiguous block row in the row */ 721 } else { 722 ierr = MatZeroRows_SeqSBAIJ_Check_Blocks(rows,is_n,bs,sizes,&bs_max);CHKERRQ(ierr); 723 } 724 ierr = ISRestoreIndices(is,&is_idx);CHKERRQ(ierr); 725 726 for (i=0,j=0; i<bs_max; j+=sizes[i],i++) { 727 row = rows[j]; /* row to be zeroed */ 728 if (row < 0 || row > A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %d out of range",row); 729 count = (sbaij->i[row/bs +1] - sbaij->i[row/bs])*bs; /* num. of elements in the row */ 730 aa = sbaij->a + sbaij->i[row/bs]*bs2 + (row%bs); 731 if (sizes[i] == bs && !sbaij->keepzeroedrows) { 732 if (diag) { 733 if (sbaij->ilen[row/bs] > 0) { 734 sbaij->ilen[row/bs] = 1; 735 sbaij->j[sbaij->i[row/bs]] = row/bs; 736 ierr = PetscMemzero(aa,count*bs*sizeof(MatScalar));CHKERRQ(ierr); 737 } 738 /* Now insert all the diagoanl values for this bs */ 739 for (k=0; k<bs; k++) { 740 ierr = (*A->ops->setvalues)(A,1,rows+j+k,1,rows+j+k,diag,INSERT_VALUES);CHKERRQ(ierr); 741 } 742 } else { /* (!diag) */ 743 sbaij->ilen[row/bs] = 0; 744 } /* end (!diag) */ 745 } else { /* (sizes[i] != bs), broken block */ 746 #if defined (PETSC_USE_DEBUG) 747 if (sizes[i] != 1) SETERRQ(1,"Internal Error. Value should be 1"); 748 #endif 749 for (k=0; k<count; k++) { 750 aa[0] = zero; 751 aa+=bs; 752 } 753 if (diag) { 754 ierr = (*A->ops->setvalues)(A,1,rows+j,1,rows+j,diag,INSERT_VALUES);CHKERRQ(ierr); 755 } 756 } 757 } 758 759 ierr = PetscFree(rows);CHKERRQ(ierr); 760 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 761 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 762 PetscFunctionReturn(0); 763 } 764 765 /* Only add/insert a(i,j) with i<=j (blocks). 766 Any a(i,j) with i>j input by user is ingored. 767 */ 768 769 #undef __FUNCT__ 770 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 771 int MatSetValues_SeqSBAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *v,InsertMode is) 772 { 773 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 774 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted=a->sorted; 775 int *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 776 int *aj=a->j,nonew=a->nonew,bs=a->bs,brow,bcol; 777 int ridx,cidx,bs2=a->bs2,ierr; 778 MatScalar *ap,value,*aa=a->a,*bap; 779 780 PetscFunctionBegin; 781 782 for (k=0; k<m; k++) { /* loop over added rows */ 783 row = im[k]; /* row number */ 784 brow = row/bs; /* block row number */ 785 if (row < 0) continue; 786 #if defined(PETSC_USE_BOPT_g) 787 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m); 788 #endif 789 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 790 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 791 rmax = imax[brow]; /* maximum space allocated for this row */ 792 nrow = ailen[brow]; /* actual length of this row */ 793 low = 0; 794 795 for (l=0; l<n; l++) { /* loop over added columns */ 796 if (in[l] < 0) continue; 797 #if defined(PETSC_USE_BOPT_g) 798 if (in[l] >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->m); 799 #endif 800 col = in[l]; 801 bcol = col/bs; /* block col number */ 802 803 if (brow <= bcol){ 804 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 805 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 806 /* element value a(k,l) */ 807 if (roworiented) { 808 value = v[l + k*n]; 809 } else { 810 value = v[k + l*m]; 811 } 812 813 /* move pointer bap to a(k,l) quickly and add/insert value */ 814 if (!sorted) low = 0; high = nrow; 815 while (high-low > 7) { 816 t = (low+high)/2; 817 if (rp[t] > bcol) high = t; 818 else low = t; 819 } 820 for (i=low; i<high; i++) { 821 /* printf("The loop of i=low.., rp[%d]=%d\n",i,rp[i]); */ 822 if (rp[i] > bcol) break; 823 if (rp[i] == bcol) { 824 bap = ap + bs2*i + bs*cidx + ridx; 825 if (is == ADD_VALUES) *bap += value; 826 else *bap = value; 827 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 828 if (brow == bcol && ridx < cidx){ 829 bap = ap + bs2*i + bs*ridx + cidx; 830 if (is == ADD_VALUES) *bap += value; 831 else *bap = value; 832 } 833 goto noinsert1; 834 } 835 } 836 837 if (nonew == 1) goto noinsert1; 838 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 839 if (nrow >= rmax) { 840 /* there is no extra room in row, therefore enlarge */ 841 int new_nz = ai[a->mbs] + CHUNKSIZE,len,*new_i,*new_j; 842 MatScalar *new_a; 843 844 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 845 846 /* Malloc new storage space */ 847 len = new_nz*(sizeof(int)+bs2*sizeof(MatScalar))+(a->mbs+1)*sizeof(int); 848 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 849 new_j = (int*)(new_a + bs2*new_nz); 850 new_i = new_j + new_nz; 851 852 /* copy over old data into new slots */ 853 for (ii=0; ii<brow+1; ii++) {new_i[ii] = ai[ii];} 854 for (ii=brow+1; ii<a->mbs+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 855 ierr = PetscMemcpy(new_j,aj,(ai[brow]+nrow)*sizeof(int));CHKERRQ(ierr); 856 len = (new_nz - CHUNKSIZE - ai[brow] - nrow); 857 ierr = PetscMemcpy(new_j+ai[brow]+nrow+CHUNKSIZE,aj+ai[brow]+nrow,len*sizeof(int));CHKERRQ(ierr); 858 ierr = PetscMemcpy(new_a,aa,(ai[brow]+nrow)*bs2*sizeof(MatScalar));CHKERRQ(ierr); 859 ierr = PetscMemzero(new_a+bs2*(ai[brow]+nrow),bs2*CHUNKSIZE*sizeof(MatScalar));CHKERRQ(ierr); 860 ierr = PetscMemcpy(new_a+bs2*(ai[brow]+nrow+CHUNKSIZE),aa+bs2*(ai[brow]+nrow),bs2*len*sizeof(MatScalar));CHKERRQ(ierr); 861 /* free up old matrix storage */ 862 ierr = PetscFree(a->a);CHKERRQ(ierr); 863 if (!a->singlemalloc) { 864 ierr = PetscFree(a->i);CHKERRQ(ierr); 865 ierr = PetscFree(a->j);CHKERRQ(ierr); 866 } 867 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 868 a->singlemalloc = PETSC_TRUE; 869 870 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 871 rmax = imax[brow] = imax[brow] + CHUNKSIZE; 872 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + bs2*sizeof(MatScalar))); 873 a->s_maxnz += bs2*CHUNKSIZE; 874 a->reallocs++; 875 a->s_nz++; 876 } 877 878 N = nrow++ - 1; 879 /* shift up all the later entries in this row */ 880 for (ii=N; ii>=i; ii--) { 881 rp[ii+1] = rp[ii]; 882 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 883 } 884 if (N>=i) { 885 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 886 } 887 rp[i] = bcol; 888 ap[bs2*i + bs*cidx + ridx] = value; 889 noinsert1:; 890 low = i; 891 /* } */ 892 } 893 } /* end of if .. if.. */ 894 } /* end of loop over added columns */ 895 ailen[brow] = nrow; 896 } /* end of loop over added rows */ 897 898 PetscFunctionReturn(0); 899 } 900 901 extern int MatCholeskyFactorSymbolic_SeqSBAIJ(Mat,IS,PetscReal,Mat*); 902 extern int MatCholeskyFactor_SeqSBAIJ(Mat,IS,PetscReal); 903 extern int MatIncreaseOverlap_SeqSBAIJ(Mat,int,IS*,int); 904 extern int MatGetSubMatrix_SeqSBAIJ(Mat,IS,IS,int,MatReuse,Mat*); 905 extern int MatGetSubMatrices_SeqSBAIJ(Mat,int,IS*,IS*,MatReuse,Mat**); 906 extern int MatMultTranspose_SeqSBAIJ(Mat,Vec,Vec); 907 extern int MatMultTransposeAdd_SeqSBAIJ(Mat,Vec,Vec,Vec); 908 extern int MatScale_SeqSBAIJ(PetscScalar*,Mat); 909 extern int MatNorm_SeqSBAIJ(Mat,NormType,PetscReal *); 910 extern int MatEqual_SeqSBAIJ(Mat,Mat,PetscTruth*); 911 extern int MatGetDiagonal_SeqSBAIJ(Mat,Vec); 912 extern int MatDiagonalScale_SeqSBAIJ(Mat,Vec,Vec); 913 extern int MatGetInfo_SeqSBAIJ(Mat,MatInfoType,MatInfo *); 914 extern int MatZeroEntries_SeqSBAIJ(Mat); 915 extern int MatGetRowMax_SeqSBAIJ(Mat,Vec); 916 917 extern int MatSolve_SeqSBAIJ_N(Mat,Vec,Vec); 918 extern int MatSolve_SeqSBAIJ_1(Mat,Vec,Vec); 919 extern int MatSolve_SeqSBAIJ_2(Mat,Vec,Vec); 920 extern int MatSolve_SeqSBAIJ_3(Mat,Vec,Vec); 921 extern int MatSolve_SeqSBAIJ_4(Mat,Vec,Vec); 922 extern int MatSolve_SeqSBAIJ_5(Mat,Vec,Vec); 923 extern int MatSolve_SeqSBAIJ_6(Mat,Vec,Vec); 924 extern int MatSolve_SeqSBAIJ_7(Mat,Vec,Vec); 925 extern int MatSolveTranspose_SeqSBAIJ_7(Mat,Vec,Vec); 926 extern int MatSolveTranspose_SeqSBAIJ_6(Mat,Vec,Vec); 927 extern int MatSolveTranspose_SeqSBAIJ_5(Mat,Vec,Vec); 928 extern int MatSolveTranspose_SeqSBAIJ_4(Mat,Vec,Vec); 929 extern int MatSolveTranspose_SeqSBAIJ_3(Mat,Vec,Vec); 930 extern int MatSolveTranspose_SeqSBAIJ_2(Mat,Vec,Vec); 931 extern int MatSolveTranspose_SeqSBAIJ_1(Mat,Vec,Vec); 932 933 extern int MatSolve_SeqSBAIJ_1_NaturalOrdering(Mat,Vec,Vec); 934 extern int MatSolve_SeqSBAIJ_2_NaturalOrdering(Mat,Vec,Vec); 935 extern int MatSolve_SeqSBAIJ_3_NaturalOrdering(Mat,Vec,Vec); 936 extern int MatSolve_SeqSBAIJ_4_NaturalOrdering(Mat,Vec,Vec); 937 extern int MatSolve_SeqSBAIJ_5_NaturalOrdering(Mat,Vec,Vec); 938 extern int MatSolve_SeqSBAIJ_6_NaturalOrdering(Mat,Vec,Vec); 939 extern int MatSolve_SeqSBAIJ_7_NaturalOrdering(Mat,Vec,Vec); 940 extern int MatSolve_SeqSBAIJ_N_NaturalOrdering(Mat,Vec,Vec); 941 942 extern int MatCholeskyFactorNumeric_SeqSBAIJ_N(Mat,Mat*); 943 extern int MatCholeskyFactorNumeric_SeqSBAIJ_1(Mat,Mat*); 944 extern int MatCholeskyFactorNumeric_SeqSBAIJ_2(Mat,Mat*); 945 extern int MatCholeskyFactorNumeric_SeqSBAIJ_3(Mat,Mat*); 946 extern int MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat,Mat*); 947 extern int MatCholeskyFactorNumeric_SeqSBAIJ_5(Mat,Mat*); 948 extern int MatCholeskyFactorNumeric_SeqSBAIJ_6(Mat,Mat*); 949 extern int MatCholeskyFactorNumeric_SeqSBAIJ_7(Mat,Mat*); 950 951 extern int MatMult_SeqSBAIJ_1(Mat,Vec,Vec); 952 extern int MatMult_SeqSBAIJ_2(Mat,Vec,Vec); 953 extern int MatMult_SeqSBAIJ_3(Mat,Vec,Vec); 954 extern int MatMult_SeqSBAIJ_4(Mat,Vec,Vec); 955 extern int MatMult_SeqSBAIJ_5(Mat,Vec,Vec); 956 extern int MatMult_SeqSBAIJ_6(Mat,Vec,Vec); 957 extern int MatMult_SeqSBAIJ_7(Mat,Vec,Vec); 958 extern int MatMult_SeqSBAIJ_N(Mat,Vec,Vec); 959 960 extern int MatMultAdd_SeqSBAIJ_1(Mat,Vec,Vec,Vec); 961 extern int MatMultAdd_SeqSBAIJ_2(Mat,Vec,Vec,Vec); 962 extern int MatMultAdd_SeqSBAIJ_3(Mat,Vec,Vec,Vec); 963 extern int MatMultAdd_SeqSBAIJ_4(Mat,Vec,Vec,Vec); 964 extern int MatMultAdd_SeqSBAIJ_5(Mat,Vec,Vec,Vec); 965 extern int MatMultAdd_SeqSBAIJ_6(Mat,Vec,Vec,Vec); 966 extern int MatMultAdd_SeqSBAIJ_7(Mat,Vec,Vec,Vec); 967 extern int MatMultAdd_SeqSBAIJ_N(Mat,Vec,Vec,Vec); 968 969 #ifdef HAVE_MatICCFactor 970 /* modefied from MatILUFactor_SeqSBAIJ, needs further work! */ 971 #undef __FUNCT__ 972 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 973 int MatICCFactor_SeqSBAIJ(Mat inA,IS row,PetscReal fill,int level) 974 { 975 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 976 Mat outA; 977 int ierr; 978 PetscTruth row_identity,col_identity; 979 980 PetscFunctionBegin; 981 /* 982 if (level != 0) SETERRQ(PETSC_ERR_SUP,"Only levels = 0 supported for in-place ILU"); 983 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 984 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 985 if (!row_identity || !col_identity) { 986 SETERRQ(1,"Row and column permutations must be identity for in-place ICC"); 987 } 988 */ 989 990 outA = inA; 991 inA->factor = FACTOR_CHOLESKY; 992 993 if (!a->diag) { 994 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 995 } 996 /* 997 Blocksize 2, 3, 4, 5, 6 and 7 have a special faster factorization/solver 998 for ILU(0) factorization with natural ordering 999 */ 1000 switch (a->bs) { 1001 case 1: 1002 inA->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering; 1003 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering; 1004 PetscLoginfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering solvetrans BS=1\n"); 1005 case 2: 1006 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1007 inA->ops->solve = MatSolve_SeqSBAIJ_2_NaturalOrdering; 1008 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_2_NaturalOrdering; 1009 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=2\n"); 1010 break; 1011 case 3: 1012 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1013 inA->ops->solve = MatSolve_SeqSBAIJ_3_NaturalOrdering; 1014 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_3_NaturalOrdering; 1015 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=3\n"); 1016 break; 1017 case 4: 1018 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1019 inA->ops->solve = MatSolve_SeqSBAIJ_4_NaturalOrdering; 1020 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_4_NaturalOrdering; 1021 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=4\n"); 1022 break; 1023 case 5: 1024 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1025 inA->ops->solve = MatSolve_SeqSBAIJ_5_NaturalOrdering; 1026 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_5_NaturalOrdering; 1027 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=5\n"); 1028 break; 1029 case 6: 1030 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1031 inA->ops->solve = MatSolve_SeqSBAIJ_6_NaturalOrdering; 1032 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_6_NaturalOrdering; 1033 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=6\n"); 1034 break; 1035 case 7: 1036 inA->ops->lufactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1037 inA->ops->solvetranspose = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1038 inA->ops->solve = MatSolve_SeqSBAIJ_7_NaturalOrdering; 1039 PetscLogInfo(inA,"MatICCFactor_SeqSBAIJ:Using special in-place natural ordering factor and solve BS=7\n"); 1040 break; 1041 default: 1042 a->row = row; 1043 a->icol = col; 1044 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1045 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1046 1047 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1048 ierr = ISInvertPermutation(col,PETSC_DECIDE, &(a->icol));CHKERRQ(ierr); 1049 PetscLogObjectParent(inA,a->icol); 1050 1051 if (!a->solve_work) { 1052 ierr = PetscMalloc((A->m+a->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1053 PetscLogObjectMemory(inA,(A->m+a->bs)*sizeof(PetscScalar)); 1054 } 1055 } 1056 1057 ierr = MatCholeskyFactorNumeric(inA,&outA);CHKERRQ(ierr); 1058 1059 PetscFunctionReturn(0); 1060 } 1061 #endif 1062 1063 #undef __FUNCT__ 1064 #define __FUNCT__ "MatPrintHelp_SeqSBAIJ" 1065 int MatPrintHelp_SeqSBAIJ(Mat A) 1066 { 1067 static PetscTruth called = PETSC_FALSE; 1068 MPI_Comm comm = A->comm; 1069 int ierr; 1070 1071 PetscFunctionBegin; 1072 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1073 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQSBAIJ and MATMPISBAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1074 ierr = (*PetscHelpPrintf)(comm," -mat_block_size <block_size>\n");CHKERRQ(ierr); 1075 PetscFunctionReturn(0); 1076 } 1077 1078 EXTERN_C_BEGIN 1079 #undef __FUNCT__ 1080 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1081 int MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,int *indices) 1082 { 1083 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1084 int i,nz,n; 1085 1086 PetscFunctionBegin; 1087 nz = baij->s_maxnz; 1088 n = mat->n; 1089 for (i=0; i<nz; i++) { 1090 baij->j[i] = indices[i]; 1091 } 1092 baij->s_nz = nz; 1093 for (i=0; i<n; i++) { 1094 baij->ilen[i] = baij->imax[i]; 1095 } 1096 1097 PetscFunctionReturn(0); 1098 } 1099 EXTERN_C_END 1100 1101 #undef __FUNCT__ 1102 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1103 /*@ 1104 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1105 in the matrix. 1106 1107 Input Parameters: 1108 + mat - the SeqSBAIJ matrix 1109 - indices - the column indices 1110 1111 Level: advanced 1112 1113 Notes: 1114 This can be called if you have precomputed the nonzero structure of the 1115 matrix and want to provide it to the matrix object to improve the performance 1116 of the MatSetValues() operation. 1117 1118 You MUST have set the correct numbers of nonzeros per row in the call to 1119 MatCreateSeqSBAIJ(). 1120 1121 MUST be called before any calls to MatSetValues() 1122 1123 .seealso: MatCreateSeqSBAIJ 1124 @*/ 1125 int MatSeqSBAIJSetColumnIndices(Mat mat,int *indices) 1126 { 1127 int ierr,(*f)(Mat,int *); 1128 1129 PetscFunctionBegin; 1130 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1131 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1132 if (f) { 1133 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1134 } else { 1135 SETERRQ(1,"Wrong type of matrix to set column indices"); 1136 } 1137 PetscFunctionReturn(0); 1138 } 1139 1140 #undef __FUNCT__ 1141 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1142 int MatSetUpPreallocation_SeqSBAIJ(Mat A) 1143 { 1144 int ierr; 1145 1146 PetscFunctionBegin; 1147 ierr = MatSeqSBAIJSetPreallocation(A,1,PETSC_DEFAULT,0);CHKERRQ(ierr); 1148 PetscFunctionReturn(0); 1149 } 1150 1151 /* -------------------------------------------------------------------*/ 1152 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1153 MatGetRow_SeqSBAIJ, 1154 MatRestoreRow_SeqSBAIJ, 1155 MatMult_SeqSBAIJ_N, 1156 MatMultAdd_SeqSBAIJ_N, 1157 MatMultTranspose_SeqSBAIJ, 1158 MatMultTransposeAdd_SeqSBAIJ, 1159 MatSolve_SeqSBAIJ_N, 1160 0, 1161 0, 1162 0, 1163 0, 1164 MatCholeskyFactor_SeqSBAIJ, 1165 MatRelax_SeqSBAIJ, 1166 MatTranspose_SeqSBAIJ, 1167 MatGetInfo_SeqSBAIJ, 1168 MatEqual_SeqSBAIJ, 1169 MatGetDiagonal_SeqSBAIJ, 1170 MatDiagonalScale_SeqSBAIJ, 1171 MatNorm_SeqSBAIJ, 1172 0, 1173 MatAssemblyEnd_SeqSBAIJ, 1174 0, 1175 MatSetOption_SeqSBAIJ, 1176 MatZeroEntries_SeqSBAIJ, 1177 MatZeroRows_SeqSBAIJ, 1178 0, 1179 0, 1180 MatCholeskyFactorSymbolic_SeqSBAIJ, 1181 MatCholeskyFactorNumeric_SeqSBAIJ_N, 1182 MatSetUpPreallocation_SeqSBAIJ, 1183 0, 1184 MatICCFactorSymbolic_SeqSBAIJ, 1185 0, 1186 0, 1187 MatDuplicate_SeqSBAIJ, 1188 0, 1189 0, 1190 0, 1191 0, 1192 0, 1193 MatGetSubMatrices_SeqSBAIJ, 1194 MatIncreaseOverlap_SeqSBAIJ, 1195 MatGetValues_SeqSBAIJ, 1196 0, 1197 MatPrintHelp_SeqSBAIJ, 1198 MatScale_SeqSBAIJ, 1199 0, 1200 0, 1201 0, 1202 MatGetBlockSize_SeqSBAIJ, 1203 MatGetRowIJ_SeqSBAIJ, 1204 MatRestoreRowIJ_SeqSBAIJ, 1205 0, 1206 0, 1207 0, 1208 0, 1209 0, 1210 0, 1211 MatSetValuesBlocked_SeqSBAIJ, 1212 MatGetSubMatrix_SeqSBAIJ, 1213 0, 1214 0, 1215 MatGetPetscMaps_Petsc, 1216 0, 1217 0, 1218 0, 1219 0, 1220 0, 1221 0, 1222 MatGetRowMax_SeqSBAIJ}; 1223 1224 EXTERN_C_BEGIN 1225 #undef __FUNCT__ 1226 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1227 int MatStoreValues_SeqSBAIJ(Mat mat) 1228 { 1229 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1230 int nz = aij->i[mat->m]*aij->bs*aij->bs2; 1231 int ierr; 1232 1233 PetscFunctionBegin; 1234 if (aij->nonew != 1) { 1235 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1236 } 1237 1238 /* allocate space for values if not already there */ 1239 if (!aij->saved_values) { 1240 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1241 } 1242 1243 /* copy values over */ 1244 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1245 PetscFunctionReturn(0); 1246 } 1247 EXTERN_C_END 1248 1249 EXTERN_C_BEGIN 1250 #undef __FUNCT__ 1251 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1252 int MatRetrieveValues_SeqSBAIJ(Mat mat) 1253 { 1254 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1255 int nz = aij->i[mat->m]*aij->bs*aij->bs2,ierr; 1256 1257 PetscFunctionBegin; 1258 if (aij->nonew != 1) { 1259 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1260 } 1261 if (!aij->saved_values) { 1262 SETERRQ(1,"Must call MatStoreValues(A);first"); 1263 } 1264 1265 /* copy values over */ 1266 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1267 PetscFunctionReturn(0); 1268 } 1269 EXTERN_C_END 1270 1271 EXTERN_C_BEGIN 1272 #undef __FUNCT__ 1273 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1274 int MatCreate_SeqSBAIJ(Mat B) 1275 { 1276 Mat_SeqSBAIJ *b; 1277 int ierr,size; 1278 1279 PetscFunctionBegin; 1280 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1281 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1282 B->m = B->M = PetscMax(B->m,B->M); 1283 B->n = B->N = PetscMax(B->n,B->N); 1284 1285 ierr = PetscNew(Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1286 B->data = (void*)b; 1287 ierr = PetscMemzero(b,sizeof(Mat_SeqSBAIJ));CHKERRQ(ierr); 1288 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1289 B->ops->destroy = MatDestroy_SeqSBAIJ; 1290 B->ops->view = MatView_SeqSBAIJ; 1291 B->factor = 0; 1292 B->lupivotthreshold = 1.0; 1293 B->mapping = 0; 1294 b->row = 0; 1295 b->icol = 0; 1296 b->reallocs = 0; 1297 b->saved_values = 0; 1298 1299 ierr = PetscMapCreateMPI(B->comm,B->m,B->M,&B->rmap);CHKERRQ(ierr); 1300 ierr = PetscMapCreateMPI(B->comm,B->n,B->N,&B->cmap);CHKERRQ(ierr); 1301 1302 b->sorted = PETSC_FALSE; 1303 b->roworiented = PETSC_TRUE; 1304 b->nonew = 0; 1305 b->diag = 0; 1306 b->solve_work = 0; 1307 b->mult_work = 0; 1308 b->spptr = 0; 1309 b->keepzeroedrows = PETSC_FALSE; 1310 1311 b->inew = 0; 1312 b->jnew = 0; 1313 b->anew = 0; 1314 b->a2anew = 0; 1315 b->permute = PETSC_FALSE; 1316 1317 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1318 "MatStoreValues_SeqSBAIJ", 1319 (void*)MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1320 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1321 "MatRetrieveValues_SeqSBAIJ", 1322 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1323 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1324 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1325 (void*)MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1326 PetscFunctionReturn(0); 1327 } 1328 EXTERN_C_END 1329 1330 #undef __FUNCT__ 1331 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1332 /*@C 1333 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1334 compressed row) format. For good matrix assembly performance the 1335 user should preallocate the matrix storage by setting the parameter nz 1336 (or the array nnz). By setting these parameters accurately, performance 1337 during matrix assembly can be increased by more than a factor of 50. 1338 1339 Collective on Mat 1340 1341 Input Parameters: 1342 + A - the symmetric matrix 1343 . bs - size of block 1344 . nz - number of block nonzeros per block row (same for all rows) 1345 - nnz - array containing the number of block nonzeros in the various block rows 1346 (possibly different for each block row) or PETSC_NULL 1347 1348 Options Database Keys: 1349 . -mat_no_unroll - uses code that does not unroll the loops in the 1350 block calculations (much slower) 1351 . -mat_block_size - size of the blocks to use 1352 1353 Level: intermediate 1354 1355 Notes: 1356 The block AIJ format is compatible with standard Fortran 77 1357 storage. That is, the stored row and column indices can begin at 1358 either one (as in Fortran) or zero. See the users' manual for details. 1359 1360 Specify the preallocated storage with either nz or nnz (not both). 1361 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1362 allocation. For additional details, see the users manual chapter on 1363 matrices. 1364 1365 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1366 @*/ 1367 int MatSeqSBAIJSetPreallocation(Mat B,int bs,int nz,int *nnz) 1368 { 1369 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1370 int i,len,ierr,mbs,bs2; 1371 PetscTruth flg; 1372 int s_nz; 1373 1374 PetscFunctionBegin; 1375 B->preallocated = PETSC_TRUE; 1376 ierr = PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1377 mbs = B->m/bs; 1378 bs2 = bs*bs; 1379 1380 if (mbs*bs != B->m) { 1381 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1382 } 1383 1384 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1385 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 1386 if (nnz) { 1387 for (i=0; i<mbs; i++) { 1388 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 1389 if (nnz[i] > mbs) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %d value %d rowlength %d",i,nnz[i],mbs); 1390 } 1391 } 1392 1393 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg);CHKERRQ(ierr); 1394 if (!flg) { 1395 switch (bs) { 1396 case 1: 1397 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1398 B->ops->solve = MatSolve_SeqSBAIJ_1; 1399 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_1; 1400 B->ops->mult = MatMult_SeqSBAIJ_1; 1401 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1402 break; 1403 case 2: 1404 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1405 B->ops->solve = MatSolve_SeqSBAIJ_2; 1406 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_2; 1407 B->ops->mult = MatMult_SeqSBAIJ_2; 1408 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1409 break; 1410 case 3: 1411 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1412 B->ops->solve = MatSolve_SeqSBAIJ_3; 1413 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_3; 1414 B->ops->mult = MatMult_SeqSBAIJ_3; 1415 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1416 break; 1417 case 4: 1418 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1419 B->ops->solve = MatSolve_SeqSBAIJ_4; 1420 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_4; 1421 B->ops->mult = MatMult_SeqSBAIJ_4; 1422 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1423 break; 1424 case 5: 1425 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1426 B->ops->solve = MatSolve_SeqSBAIJ_5; 1427 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_5; 1428 B->ops->mult = MatMult_SeqSBAIJ_5; 1429 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1430 break; 1431 case 6: 1432 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1433 B->ops->solve = MatSolve_SeqSBAIJ_6; 1434 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_6; 1435 B->ops->mult = MatMult_SeqSBAIJ_6; 1436 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1437 break; 1438 case 7: 1439 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1440 B->ops->solve = MatSolve_SeqSBAIJ_7; 1441 B->ops->solvetranspose = MatSolveTranspose_SeqSBAIJ_7; 1442 B->ops->mult = MatMult_SeqSBAIJ_7; 1443 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1444 break; 1445 } 1446 } 1447 1448 b->mbs = mbs; 1449 b->nbs = mbs; 1450 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 1451 if (!nnz) { 1452 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1453 else if (nz <= 0) nz = 1; 1454 for (i=0; i<mbs; i++) { 1455 b->imax[i] = nz; 1456 } 1457 nz = nz*mbs; /* total nz */ 1458 } else { 1459 nz = 0; 1460 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1461 } 1462 /* s_nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1463 s_nz = nz; 1464 1465 /* allocate the matrix space */ 1466 len = s_nz*sizeof(int) + s_nz*bs2*sizeof(MatScalar) + (B->m+1)*sizeof(int); 1467 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 1468 ierr = PetscMemzero(b->a,s_nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1469 b->j = (int*)(b->a + s_nz*bs2); 1470 ierr = PetscMemzero(b->j,s_nz*sizeof(int));CHKERRQ(ierr); 1471 b->i = b->j + s_nz; 1472 b->singlemalloc = PETSC_TRUE; 1473 1474 /* pointer to beginning of each row */ 1475 b->i[0] = 0; 1476 for (i=1; i<mbs+1; i++) { 1477 b->i[i] = b->i[i-1] + b->imax[i-1]; 1478 } 1479 1480 /* b->ilen will count nonzeros in each block row so far. */ 1481 ierr = PetscMalloc((mbs+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 1482 PetscLogObjectMemory(B,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1483 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1484 1485 b->bs = bs; 1486 b->bs2 = bs2; 1487 b->s_nz = 0; 1488 b->s_maxnz = s_nz*bs2; 1489 1490 b->inew = 0; 1491 b->jnew = 0; 1492 b->anew = 0; 1493 b->a2anew = 0; 1494 b->permute = PETSC_FALSE; 1495 PetscFunctionReturn(0); 1496 } 1497 1498 1499 #undef __FUNCT__ 1500 #define __FUNCT__ "MatCreateSeqSBAIJ" 1501 /*@C 1502 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1503 compressed row) format. For good matrix assembly performance the 1504 user should preallocate the matrix storage by setting the parameter nz 1505 (or the array nnz). By setting these parameters accurately, performance 1506 during matrix assembly can be increased by more than a factor of 50. 1507 1508 Collective on MPI_Comm 1509 1510 Input Parameters: 1511 + comm - MPI communicator, set to PETSC_COMM_SELF 1512 . bs - size of block 1513 . m - number of rows, or number of columns 1514 . nz - number of block nonzeros per block row (same for all rows) 1515 - nnz - array containing the number of block nonzeros in the various block rows 1516 (possibly different for each block row) or PETSC_NULL 1517 1518 Output Parameter: 1519 . A - the symmetric matrix 1520 1521 Options Database Keys: 1522 . -mat_no_unroll - uses code that does not unroll the loops in the 1523 block calculations (much slower) 1524 . -mat_block_size - size of the blocks to use 1525 1526 Level: intermediate 1527 1528 Notes: 1529 The block AIJ format is compatible with standard Fortran 77 1530 storage. That is, the stored row and column indices can begin at 1531 either one (as in Fortran) or zero. See the users' manual for details. 1532 1533 Specify the preallocated storage with either nz or nnz (not both). 1534 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1535 allocation. For additional details, see the users manual chapter on 1536 matrices. 1537 1538 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1539 @*/ 1540 int MatCreateSeqSBAIJ(MPI_Comm comm,int bs,int m,int n,int nz,int *nnz,Mat *A) 1541 { 1542 int ierr; 1543 1544 PetscFunctionBegin; 1545 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1546 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1547 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,nnz);CHKERRQ(ierr); 1548 PetscFunctionReturn(0); 1549 } 1550 1551 #undef __FUNCT__ 1552 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1553 int MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1554 { 1555 Mat C; 1556 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1557 int i,len,mbs = a->mbs,nz = a->s_nz,bs2 =a->bs2,ierr; 1558 1559 PetscFunctionBegin; 1560 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1561 1562 *B = 0; 1563 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 1564 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 1565 c = (Mat_SeqSBAIJ*)C->data; 1566 1567 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1568 C->preallocated = PETSC_TRUE; 1569 C->factor = A->factor; 1570 c->row = 0; 1571 c->icol = 0; 1572 c->saved_values = 0; 1573 c->keepzeroedrows = a->keepzeroedrows; 1574 C->assembled = PETSC_TRUE; 1575 1576 c->bs = a->bs; 1577 c->bs2 = a->bs2; 1578 c->mbs = a->mbs; 1579 c->nbs = a->nbs; 1580 1581 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 1582 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 1583 for (i=0; i<mbs; i++) { 1584 c->imax[i] = a->imax[i]; 1585 c->ilen[i] = a->ilen[i]; 1586 } 1587 1588 /* allocate the matrix space */ 1589 c->singlemalloc = PETSC_TRUE; 1590 len = (mbs+1)*sizeof(int) + nz*(bs2*sizeof(MatScalar) + sizeof(int)); 1591 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 1592 c->j = (int*)(c->a + nz*bs2); 1593 c->i = c->j + nz; 1594 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(int));CHKERRQ(ierr); 1595 if (mbs > 0) { 1596 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(int));CHKERRQ(ierr); 1597 if (cpvalues == MAT_COPY_VALUES) { 1598 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1599 } else { 1600 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1601 } 1602 } 1603 1604 PetscLogObjectMemory(C,len+2*(mbs+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqSBAIJ)); 1605 c->sorted = a->sorted; 1606 c->roworiented = a->roworiented; 1607 c->nonew = a->nonew; 1608 1609 if (a->diag) { 1610 ierr = PetscMalloc((mbs+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 1611 PetscLogObjectMemory(C,(mbs+1)*sizeof(int)); 1612 for (i=0; i<mbs; i++) { 1613 c->diag[i] = a->diag[i]; 1614 } 1615 } else c->diag = 0; 1616 c->s_nz = a->s_nz; 1617 c->s_maxnz = a->s_maxnz; 1618 c->solve_work = 0; 1619 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1620 c->mult_work = 0; 1621 *B = C; 1622 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 1623 PetscFunctionReturn(0); 1624 } 1625 1626 EXTERN_C_BEGIN 1627 #undef __FUNCT__ 1628 #define __FUNCT__ "MatLoad_SeqSBAIJ" 1629 int MatLoad_SeqSBAIJ(PetscViewer viewer,MatType type,Mat *A) 1630 { 1631 Mat_SeqSBAIJ *a; 1632 Mat B; 1633 int i,nz,ierr,fd,header[4],size,*rowlengths=0,M,N,bs=1; 1634 int *mask,mbs,*jj,j,rowcount,nzcount,k,*browlengths,*s_browlengths,maskcount; 1635 int kmax,jcount,block,idx,point,nzcountb,extra_rows; 1636 int *masked,nmask,tmp,bs2,ishift; 1637 PetscScalar *aa; 1638 MPI_Comm comm = ((PetscObject)viewer)->comm; 1639 1640 PetscFunctionBegin; 1641 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1642 bs2 = bs*bs; 1643 1644 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1645 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1646 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1647 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1648 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1649 M = header[1]; N = header[2]; nz = header[3]; 1650 1651 if (header[3] < 0) { 1652 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1653 } 1654 1655 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1656 1657 /* 1658 This code adds extra rows to make sure the number of rows is 1659 divisible by the blocksize 1660 */ 1661 mbs = M/bs; 1662 extra_rows = bs - M + bs*(mbs); 1663 if (extra_rows == bs) extra_rows = 0; 1664 else mbs++; 1665 if (extra_rows) { 1666 PetscLogInfo(0,"MatLoad_SeqSBAIJ:Padding loaded matrix to match blocksize\n"); 1667 } 1668 1669 /* read in row lengths */ 1670 ierr = PetscMalloc((M+extra_rows)*sizeof(int),&rowlengths);CHKERRQ(ierr); 1671 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1672 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1673 1674 /* read in column indices */ 1675 ierr = PetscMalloc((nz+extra_rows)*sizeof(int),&jj);CHKERRQ(ierr); 1676 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1677 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1678 1679 /* loop over row lengths determining block row lengths */ 1680 ierr = PetscMalloc(mbs*sizeof(int),&browlengths);CHKERRQ(ierr); 1681 ierr = PetscMalloc(mbs*sizeof(int),&s_browlengths);CHKERRQ(ierr); 1682 ierr = PetscMemzero(s_browlengths,mbs*sizeof(int));CHKERRQ(ierr); 1683 ierr = PetscMalloc(2*mbs*sizeof(int),&mask);CHKERRQ(ierr); 1684 ierr = PetscMemzero(mask,mbs*sizeof(int));CHKERRQ(ierr); 1685 masked = mask + mbs; 1686 rowcount = 0; nzcount = 0; 1687 for (i=0; i<mbs; i++) { 1688 nmask = 0; 1689 for (j=0; j<bs; j++) { 1690 kmax = rowlengths[rowcount]; 1691 for (k=0; k<kmax; k++) { 1692 tmp = jj[nzcount++]/bs; /* block col. index */ 1693 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1694 } 1695 rowcount++; 1696 } 1697 s_browlengths[i] += nmask; 1698 browlengths[i] = 2*s_browlengths[i]; 1699 1700 /* zero out the mask elements we set */ 1701 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1702 } 1703 1704 /* create our matrix */ 1705 ierr = MatCreateSeqSBAIJ(comm,bs,M+extra_rows,N+extra_rows,0,browlengths,A);CHKERRQ(ierr); 1706 B = *A; 1707 a = (Mat_SeqSBAIJ*)B->data; 1708 1709 /* set matrix "i" values */ 1710 a->i[0] = 0; 1711 for (i=1; i<= mbs; i++) { 1712 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1713 a->ilen[i-1] = s_browlengths[i-1]; 1714 } 1715 a->s_nz = 0; 1716 for (i=0; i<mbs; i++) a->s_nz += s_browlengths[i]; 1717 1718 /* read in nonzero values */ 1719 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 1720 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1721 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1722 1723 /* set "a" and "j" values into matrix */ 1724 nzcount = 0; jcount = 0; 1725 for (i=0; i<mbs; i++) { 1726 nzcountb = nzcount; 1727 nmask = 0; 1728 for (j=0; j<bs; j++) { 1729 kmax = rowlengths[i*bs+j]; 1730 for (k=0; k<kmax; k++) { 1731 tmp = jj[nzcount++]/bs; /* block col. index */ 1732 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1733 } 1734 } 1735 /* sort the masked values */ 1736 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1737 1738 /* set "j" values into matrix */ 1739 maskcount = 1; 1740 for (j=0; j<nmask; j++) { 1741 a->j[jcount++] = masked[j]; 1742 mask[masked[j]] = maskcount++; 1743 } 1744 1745 /* set "a" values into matrix */ 1746 ishift = bs2*a->i[i]; 1747 for (j=0; j<bs; j++) { 1748 kmax = rowlengths[i*bs+j]; 1749 for (k=0; k<kmax; k++) { 1750 tmp = jj[nzcountb]/bs ; /* block col. index */ 1751 if (tmp >= i){ 1752 block = mask[tmp] - 1; 1753 point = jj[nzcountb] - bs*tmp; 1754 idx = ishift + bs2*block + j + bs*point; 1755 a->a[idx] = aa[nzcountb]; 1756 } 1757 nzcountb++; 1758 } 1759 } 1760 /* zero out the mask elements we set */ 1761 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1762 } 1763 if (jcount != a->s_nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 1764 1765 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1766 ierr = PetscFree(browlengths);CHKERRQ(ierr); 1767 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 1768 ierr = PetscFree(aa);CHKERRQ(ierr); 1769 ierr = PetscFree(jj);CHKERRQ(ierr); 1770 ierr = PetscFree(mask);CHKERRQ(ierr); 1771 1772 B->assembled = PETSC_TRUE; 1773 ierr = MatView_Private(B);CHKERRQ(ierr); 1774 PetscFunctionReturn(0); 1775 } 1776 EXTERN_C_END 1777 1778 #undef __FUNCT__ 1779 #define __FUNCT__ "MatRelax_SeqSBAIJ" 1780 int MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,int lits,Vec xx) 1781 { 1782 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1783 MatScalar *aa=a->a,*v,*v1; 1784 PetscScalar *x,*b,*t,sum,d; 1785 int m=a->mbs,bs=a->bs,*ai=a->i,*aj=a->j,ierr; 1786 int nz,nz1,*vj,*vj1,i; 1787 1788 PetscFunctionBegin; 1789 its = its*lits; 1790 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits); 1791 1792 if (bs > 1) 1793 SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented"); 1794 1795 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1796 if (xx != bb) { 1797 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 1798 } else { 1799 b = x; 1800 } 1801 1802 ierr = PetscMalloc(m*sizeof(PetscScalar),&t);CHKERRQ(ierr); 1803 1804 if (flag & SOR_ZERO_INITIAL_GUESS) { 1805 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1806 for (i=0; i<m; i++) 1807 t[i] = b[i]; 1808 1809 for (i=0; i<m; i++){ 1810 d = *(aa + ai[i]); /* diag[i] */ 1811 v = aa + ai[i] + 1; 1812 vj = aj + ai[i] + 1; 1813 nz = ai[i+1] - ai[i] - 1; 1814 x[i] = omega*t[i]/d; 1815 while (nz--) t[*vj++] -= x[i]*(*v++); /* update rhs */ 1816 PetscLogFlops(2*nz-1); 1817 } 1818 } 1819 1820 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1821 for (i=0; i<m; i++) 1822 t[i] = b[i]; 1823 1824 for (i=0; i<m-1; i++){ /* update rhs */ 1825 v = aa + ai[i] + 1; 1826 vj = aj + ai[i] + 1; 1827 nz = ai[i+1] - ai[i] - 1; 1828 while (nz--) t[*vj++] -= x[i]*(*v++); 1829 PetscLogFlops(2*nz-1); 1830 } 1831 for (i=m-1; i>=0; i--){ 1832 d = *(aa + ai[i]); 1833 v = aa + ai[i] + 1; 1834 vj = aj + ai[i] + 1; 1835 nz = ai[i+1] - ai[i] - 1; 1836 sum = t[i]; 1837 while (nz--) sum -= x[*vj++]*(*v++); 1838 PetscLogFlops(2*nz-1); 1839 x[i] = (1-omega)*x[i] + omega*sum/d; 1840 } 1841 } 1842 its--; 1843 } 1844 1845 while (its--) { 1846 /* 1847 forward sweep: 1848 for i=0,...,m-1: 1849 sum[i] = (b[i] - U(i,:)x )/d[i]; 1850 x[i] = (1-omega)x[i] + omega*sum[i]; 1851 b = b - x[i]*U^T(i,:); 1852 1853 */ 1854 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1855 for (i=0; i<m; i++) 1856 t[i] = b[i]; 1857 1858 for (i=0; i<m; i++){ 1859 d = *(aa + ai[i]); /* diag[i] */ 1860 v = aa + ai[i] + 1; v1=v; 1861 vj = aj + ai[i] + 1; vj1=vj; 1862 nz = ai[i+1] - ai[i] - 1; nz1=nz; 1863 sum = t[i]; 1864 while (nz1--) sum -= (*v1++)*x[*vj1++]; 1865 x[i] = (1-omega)*x[i] + omega*sum/d; 1866 while (nz--) t[*vj++] -= x[i]*(*v++); 1867 PetscLogFlops(4*nz-2); 1868 } 1869 } 1870 1871 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1872 /* 1873 backward sweep: 1874 b = b - x[i]*U^T(i,:), i=0,...,n-2 1875 for i=m-1,...,0: 1876 sum[i] = (b[i] - U(i,:)x )/d[i]; 1877 x[i] = (1-omega)x[i] + omega*sum[i]; 1878 */ 1879 for (i=0; i<m; i++) 1880 t[i] = b[i]; 1881 1882 for (i=0; i<m-1; i++){ /* update rhs */ 1883 v = aa + ai[i] + 1; 1884 vj = aj + ai[i] + 1; 1885 nz = ai[i+1] - ai[i] - 1; 1886 while (nz--) t[*vj++] -= x[i]*(*v++); 1887 PetscLogFlops(2*nz-1); 1888 } 1889 for (i=m-1; i>=0; i--){ 1890 d = *(aa + ai[i]); 1891 v = aa + ai[i] + 1; 1892 vj = aj + ai[i] + 1; 1893 nz = ai[i+1] - ai[i] - 1; 1894 sum = t[i]; 1895 while (nz--) sum -= x[*vj++]*(*v++); 1896 PetscLogFlops(2*nz-1); 1897 x[i] = (1-omega)*x[i] + omega*sum/d; 1898 } 1899 } 1900 } 1901 1902 ierr = PetscFree(t); CHKERRQ(ierr); 1903 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1904 if (bb != xx) { 1905 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 1906 } 1907 1908 PetscFunctionReturn(0); 1909 } 1910 1911 1912 1913 1914 1915 1916