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