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