1 2 /* 3 Defines the basic matrix operations for the SBAIJ (compressed row) 4 matrix storage format. 5 */ 6 #include <../src/mat/impls/baij/seq/baij.h> /*I "petscmat.h" I*/ 7 #include <../src/mat/impls/sbaij/seq/sbaij.h> 8 #include <petscblaslapack.h> 9 10 #include <../src/mat/impls/sbaij/seq/relax.h> 11 #define USESHORT 12 #include <../src/mat/impls/sbaij/seq/relax.h> 13 14 #if defined(PETSC_HAVE_ELEMENTAL) 15 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 16 #endif 17 #if defined(PETSC_HAVE_SCALAPACK) 18 PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*); 19 #endif 20 PETSC_INTERN PetscErrorCode MatConvert_MPISBAIJ_Basic(Mat,MatType,MatReuse,Mat*); 21 22 /* 23 Checks for missing diagonals 24 */ 25 PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A,PetscBool *missing,PetscInt *dd) 26 { 27 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 28 PetscErrorCode ierr; 29 PetscInt *diag,*ii = a->i,i; 30 31 PetscFunctionBegin; 32 ierr = MatMarkDiagonal_SeqSBAIJ(A);CHKERRQ(ierr); 33 *missing = PETSC_FALSE; 34 if (A->rmap->n > 0 && !ii) { 35 *missing = PETSC_TRUE; 36 if (dd) *dd = 0; 37 ierr = PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n");CHKERRQ(ierr); 38 } else { 39 diag = a->diag; 40 for (i=0; i<a->mbs; i++) { 41 if (diag[i] >= ii[i+1]) { 42 *missing = PETSC_TRUE; 43 if (dd) *dd = i; 44 break; 45 } 46 } 47 } 48 PetscFunctionReturn(0); 49 } 50 51 PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A) 52 { 53 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 54 PetscErrorCode ierr; 55 PetscInt i,j; 56 57 PetscFunctionBegin; 58 if (!a->diag) { 59 ierr = PetscMalloc1(a->mbs,&a->diag);CHKERRQ(ierr); 60 ierr = PetscLogObjectMemory((PetscObject)A,a->mbs*sizeof(PetscInt));CHKERRQ(ierr); 61 a->free_diag = PETSC_TRUE; 62 } 63 for (i=0; i<a->mbs; i++) { 64 a->diag[i] = a->i[i+1]; 65 for (j=a->i[i]; j<a->i[i+1]; j++) { 66 if (a->j[j] == i) { 67 a->diag[i] = j; 68 break; 69 } 70 } 71 } 72 PetscFunctionReturn(0); 73 } 74 75 static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *inia[],const PetscInt *inja[],PetscBool *done) 76 { 77 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 78 PetscErrorCode ierr; 79 PetscInt i,j,n = a->mbs,nz = a->i[n],*tia,*tja,bs = A->rmap->bs,k,l,cnt; 80 PetscInt **ia = (PetscInt**)inia,**ja = (PetscInt**)inja; 81 82 PetscFunctionBegin; 83 *nn = n; 84 if (!ia) PetscFunctionReturn(0); 85 if (symmetric) { 86 ierr = MatToSymmetricIJ_SeqAIJ(n,a->i,a->j,PETSC_FALSE,0,0,&tia,&tja);CHKERRQ(ierr); 87 nz = tia[n]; 88 } else { 89 tia = a->i; tja = a->j; 90 } 91 92 if (!blockcompressed && bs > 1) { 93 (*nn) *= bs; 94 /* malloc & create the natural set of indices */ 95 ierr = PetscMalloc1((n+1)*bs,ia);CHKERRQ(ierr); 96 if (n) { 97 (*ia)[0] = oshift; 98 for (j=1; j<bs; j++) { 99 (*ia)[j] = (tia[1]-tia[0])*bs+(*ia)[j-1]; 100 } 101 } 102 103 for (i=1; i<n; i++) { 104 (*ia)[i*bs] = (tia[i]-tia[i-1])*bs + (*ia)[i*bs-1]; 105 for (j=1; j<bs; j++) { 106 (*ia)[i*bs+j] = (tia[i+1]-tia[i])*bs + (*ia)[i*bs+j-1]; 107 } 108 } 109 if (n) { 110 (*ia)[n*bs] = (tia[n]-tia[n-1])*bs + (*ia)[n*bs-1]; 111 } 112 113 if (inja) { 114 ierr = PetscMalloc1(nz*bs*bs,ja);CHKERRQ(ierr); 115 cnt = 0; 116 for (i=0; i<n; i++) { 117 for (j=0; j<bs; j++) { 118 for (k=tia[i]; k<tia[i+1]; k++) { 119 for (l=0; l<bs; l++) { 120 (*ja)[cnt++] = bs*tja[k] + l; 121 } 122 } 123 } 124 } 125 } 126 127 if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */ 128 ierr = PetscFree(tia);CHKERRQ(ierr); 129 ierr = PetscFree(tja);CHKERRQ(ierr); 130 } 131 } else if (oshift == 1) { 132 if (symmetric) { 133 nz = tia[A->rmap->n/bs]; 134 /* add 1 to i and j indices */ 135 for (i=0; i<A->rmap->n/bs+1; i++) tia[i] = tia[i] + 1; 136 *ia = tia; 137 if (ja) { 138 for (i=0; i<nz; i++) tja[i] = tja[i] + 1; 139 *ja = tja; 140 } 141 } else { 142 nz = a->i[A->rmap->n/bs]; 143 /* malloc space and add 1 to i and j indices */ 144 ierr = PetscMalloc1(A->rmap->n/bs+1,ia);CHKERRQ(ierr); 145 for (i=0; i<A->rmap->n/bs+1; i++) (*ia)[i] = a->i[i] + 1; 146 if (ja) { 147 ierr = PetscMalloc1(nz,ja);CHKERRQ(ierr); 148 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 149 } 150 } 151 } else { 152 *ia = tia; 153 if (ja) *ja = tja; 154 } 155 PetscFunctionReturn(0); 156 } 157 158 static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool blockcompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 159 { 160 PetscErrorCode ierr; 161 162 PetscFunctionBegin; 163 if (!ia) PetscFunctionReturn(0); 164 if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) { 165 ierr = PetscFree(*ia);CHKERRQ(ierr); 166 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 167 } 168 PetscFunctionReturn(0); 169 } 170 171 PetscErrorCode MatDestroy_SeqSBAIJ(Mat A) 172 { 173 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 #if defined(PETSC_USE_LOG) 178 PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap->N,a->nz); 179 #endif 180 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 181 if (a->free_diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 182 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 183 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 184 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 185 ierr = PetscFree(a->idiag);CHKERRQ(ierr); 186 ierr = PetscFree(a->inode.size);CHKERRQ(ierr); 187 if (a->free_imax_ilen) {ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr);} 188 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 189 ierr = PetscFree(a->sor_work);CHKERRQ(ierr); 190 ierr = PetscFree(a->solves_work);CHKERRQ(ierr); 191 ierr = PetscFree(a->mult_work);CHKERRQ(ierr); 192 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 193 if (a->free_jshort) {ierr = PetscFree(a->jshort);CHKERRQ(ierr);} 194 ierr = PetscFree(a->inew);CHKERRQ(ierr); 195 ierr = MatDestroy(&a->parent);CHKERRQ(ierr); 196 ierr = PetscFree(A->data);CHKERRQ(ierr); 197 198 ierr = PetscObjectChangeTypeName((PetscObject)A,NULL);CHKERRQ(ierr); 199 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);CHKERRQ(ierr); 200 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);CHKERRQ(ierr); 201 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C",NULL);CHKERRQ(ierr); 202 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C",NULL);CHKERRQ(ierr); 203 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C",NULL);CHKERRQ(ierr); 204 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",NULL);CHKERRQ(ierr); 205 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr); 206 #if defined(PETSC_HAVE_ELEMENTAL) 207 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_elemental_C",NULL);CHKERRQ(ierr); 208 #endif 209 #if defined(PETSC_HAVE_SCALAPACK) 210 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_scalapack_C",NULL);CHKERRQ(ierr); 211 #endif 212 PetscFunctionReturn(0); 213 } 214 215 PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscBool flg) 216 { 217 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 218 #if defined(PETSC_USE_COMPLEX) 219 PetscInt bs; 220 #endif 221 PetscErrorCode ierr; 222 223 PetscFunctionBegin; 224 #if defined(PETSC_USE_COMPLEX) 225 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 226 #endif 227 switch (op) { 228 case MAT_ROW_ORIENTED: 229 a->roworiented = flg; 230 break; 231 case MAT_KEEP_NONZERO_PATTERN: 232 a->keepnonzeropattern = flg; 233 break; 234 case MAT_NEW_NONZERO_LOCATIONS: 235 a->nonew = (flg ? 0 : 1); 236 break; 237 case MAT_NEW_NONZERO_LOCATION_ERR: 238 a->nonew = (flg ? -1 : 0); 239 break; 240 case MAT_NEW_NONZERO_ALLOCATION_ERR: 241 a->nonew = (flg ? -2 : 0); 242 break; 243 case MAT_UNUSED_NONZERO_LOCATION_ERR: 244 a->nounused = (flg ? -1 : 0); 245 break; 246 case MAT_FORCE_DIAGONAL_ENTRIES: 247 case MAT_IGNORE_OFF_PROC_ENTRIES: 248 case MAT_USE_HASH_TABLE: 249 case MAT_SORTED_FULL: 250 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 251 break; 252 case MAT_HERMITIAN: 253 #if defined(PETSC_USE_COMPLEX) 254 if (flg) { /* disable transpose ops */ 255 if (bs > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for Hermitian with block size greater than 1"); 256 A->ops->multtranspose = NULL; 257 A->ops->multtransposeadd = NULL; 258 A->symmetric = PETSC_FALSE; 259 } 260 #endif 261 break; 262 case MAT_SYMMETRIC: 263 case MAT_SPD: 264 #if defined(PETSC_USE_COMPLEX) 265 if (flg) { /* An hermitian and symmetric matrix has zero imaginary part (restore back transpose ops) */ 266 A->ops->multtranspose = A->ops->mult; 267 A->ops->multtransposeadd = A->ops->multadd; 268 } 269 #endif 270 break; 271 /* These options are handled directly by MatSetOption() */ 272 case MAT_STRUCTURALLY_SYMMETRIC: 273 case MAT_SYMMETRY_ETERNAL: 274 case MAT_STRUCTURE_ONLY: 275 /* These options are handled directly by MatSetOption() */ 276 break; 277 case MAT_IGNORE_LOWER_TRIANGULAR: 278 a->ignore_ltriangular = flg; 279 break; 280 case MAT_ERROR_LOWER_TRIANGULAR: 281 a->ignore_ltriangular = flg; 282 break; 283 case MAT_GETROW_UPPERTRIANGULAR: 284 a->getrow_utriangular = flg; 285 break; 286 case MAT_SUBMAT_SINGLEIS: 287 break; 288 default: 289 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 290 } 291 PetscFunctionReturn(0); 292 } 293 294 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 295 { 296 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 if (A && !a->getrow_utriangular) SETERRQ(PETSC_COMM_SELF,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,PETSC_TRUE) or MatGetRowUpperTriangular()"); 301 302 /* Get the upper triangular part of the row */ 303 ierr = MatGetRow_SeqBAIJ_private(A,row,nz,idx,v,a->i,a->j,a->a);CHKERRQ(ierr); 304 PetscFunctionReturn(0); 305 } 306 307 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 308 { 309 PetscErrorCode ierr; 310 311 PetscFunctionBegin; 312 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);} 313 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);} 314 PetscFunctionReturn(0); 315 } 316 317 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 318 { 319 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 320 321 PetscFunctionBegin; 322 a->getrow_utriangular = PETSC_TRUE; 323 PetscFunctionReturn(0); 324 } 325 326 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) 327 { 328 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 329 330 PetscFunctionBegin; 331 a->getrow_utriangular = PETSC_FALSE; 332 PetscFunctionReturn(0); 333 } 334 335 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B) 336 { 337 PetscErrorCode ierr; 338 339 PetscFunctionBegin; 340 if (reuse == MAT_INITIAL_MATRIX) { 341 ierr = MatDuplicate(A,MAT_COPY_VALUES,B);CHKERRQ(ierr); 342 } else if (reuse == MAT_REUSE_MATRIX) { 343 ierr = MatCopy(A,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 344 } 345 PetscFunctionReturn(0); 346 } 347 348 PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 349 { 350 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 351 PetscErrorCode ierr; 352 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2; 353 PetscViewerFormat format; 354 PetscInt *diag; 355 356 PetscFunctionBegin; 357 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 358 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 359 ierr = PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);CHKERRQ(ierr); 360 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 361 Mat aij; 362 const char *matname; 363 364 if (A->factortype && bs>1) { 365 ierr = PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");CHKERRQ(ierr); 366 PetscFunctionReturn(0); 367 } 368 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);CHKERRQ(ierr); 369 ierr = PetscObjectGetName((PetscObject)A,&matname);CHKERRQ(ierr); 370 ierr = PetscObjectSetName((PetscObject)aij,matname);CHKERRQ(ierr); 371 ierr = MatView(aij,viewer);CHKERRQ(ierr); 372 ierr = MatDestroy(&aij);CHKERRQ(ierr); 373 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 374 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 375 for (i=0; i<a->mbs; i++) { 376 for (j=0; j<bs; j++) { 377 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 378 for (k=a->i[i]; k<a->i[i+1]; k++) { 379 for (l=0; l<bs; l++) { 380 #if defined(PETSC_USE_COMPLEX) 381 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 382 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l, 383 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 384 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 385 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l, 386 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 387 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 388 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 389 } 390 #else 391 if (a->a[bs2*k + l*bs + j] != 0.0) { 392 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 393 } 394 #endif 395 } 396 } 397 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 398 } 399 } 400 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 401 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 402 PetscFunctionReturn(0); 403 } else { 404 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 405 if (A->factortype) { /* for factored matrix */ 406 if (bs>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet"); 407 408 diag=a->diag; 409 for (i=0; i<a->mbs; i++) { /* for row block i */ 410 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 411 /* diagonal entry */ 412 #if defined(PETSC_USE_COMPLEX) 413 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) { 414 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),(double)PetscImaginaryPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 415 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) { 416 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]),-(double)PetscImaginaryPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 417 } else { 418 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 419 } 420 #else 421 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]));CHKERRQ(ierr); 422 #endif 423 /* off-diagonal entries */ 424 for (k=a->i[i]; k<a->i[i+1]-1; k++) { 425 #if defined(PETSC_USE_COMPLEX) 426 if (PetscImaginaryPart(a->a[k]) > 0.0) { 427 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 428 } else if (PetscImaginaryPart(a->a[k]) < 0.0) { 429 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 430 } else { 431 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]));CHKERRQ(ierr); 432 } 433 #else 434 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[k],(double)a->a[k]);CHKERRQ(ierr); 435 #endif 436 } 437 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 438 } 439 440 } else { /* for non-factored matrix */ 441 for (i=0; i<a->mbs; i++) { /* for row block i */ 442 for (j=0; j<bs; j++) { /* for row bs*i + j */ 443 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 444 for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */ 445 for (l=0; l<bs; l++) { /* for column */ 446 #if defined(PETSC_USE_COMPLEX) 447 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 448 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l, 449 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 450 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 451 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l, 452 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 453 } else { 454 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 455 } 456 #else 457 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 458 #endif 459 } 460 } 461 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 462 } 463 } 464 } 465 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 466 } 467 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 468 PetscFunctionReturn(0); 469 } 470 471 #include <petscdraw.h> 472 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 473 { 474 Mat A = (Mat) Aa; 475 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 476 PetscErrorCode ierr; 477 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 478 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 479 MatScalar *aa; 480 PetscViewer viewer; 481 482 PetscFunctionBegin; 483 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 484 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 485 486 /* loop over matrix elements drawing boxes */ 487 488 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 489 ierr = PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric");CHKERRQ(ierr); 490 /* Blue for negative, Cyan for zero and Red for positive */ 491 color = PETSC_DRAW_BLUE; 492 for (i=0,row=0; i<mbs; i++,row+=bs) { 493 for (j=a->i[i]; j<a->i[i+1]; j++) { 494 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 495 x_l = a->j[j]*bs; x_r = x_l + 1.0; 496 aa = a->a + j*bs2; 497 for (k=0; k<bs; k++) { 498 for (l=0; l<bs; l++) { 499 if (PetscRealPart(*aa++) >= 0.) continue; 500 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 501 } 502 } 503 } 504 } 505 color = PETSC_DRAW_CYAN; 506 for (i=0,row=0; i<mbs; i++,row+=bs) { 507 for (j=a->i[i]; j<a->i[i+1]; j++) { 508 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 509 x_l = a->j[j]*bs; x_r = x_l + 1.0; 510 aa = a->a + j*bs2; 511 for (k=0; k<bs; k++) { 512 for (l=0; l<bs; l++) { 513 if (PetscRealPart(*aa++) != 0.) continue; 514 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 515 } 516 } 517 } 518 } 519 color = PETSC_DRAW_RED; 520 for (i=0,row=0; i<mbs; i++,row+=bs) { 521 for (j=a->i[i]; j<a->i[i+1]; j++) { 522 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 523 x_l = a->j[j]*bs; x_r = x_l + 1.0; 524 aa = a->a + j*bs2; 525 for (k=0; k<bs; k++) { 526 for (l=0; l<bs; l++) { 527 if (PetscRealPart(*aa++) <= 0.) continue; 528 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 529 } 530 } 531 } 532 } 533 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 534 PetscFunctionReturn(0); 535 } 536 537 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 538 { 539 PetscErrorCode ierr; 540 PetscReal xl,yl,xr,yr,w,h; 541 PetscDraw draw; 542 PetscBool isnull; 543 544 PetscFunctionBegin; 545 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 546 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 547 if (isnull) PetscFunctionReturn(0); 548 549 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 550 xr += w; yr += h; xl = -w; yl = -h; 551 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 552 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 553 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 554 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 555 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 556 PetscFunctionReturn(0); 557 } 558 559 /* Used for both MPIBAIJ and MPISBAIJ matrices */ 560 #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary 561 562 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 563 { 564 PetscErrorCode ierr; 565 PetscBool iascii,isbinary,isdraw; 566 567 PetscFunctionBegin; 568 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 569 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 570 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 571 if (iascii) { 572 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 573 } else if (isbinary) { 574 ierr = MatView_SeqSBAIJ_Binary(A,viewer);CHKERRQ(ierr); 575 } else if (isdraw) { 576 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 577 } else { 578 Mat B; 579 const char *matname; 580 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 581 ierr = PetscObjectGetName((PetscObject)A,&matname);CHKERRQ(ierr); 582 ierr = PetscObjectSetName((PetscObject)B,matname);CHKERRQ(ierr); 583 ierr = MatView(B,viewer);CHKERRQ(ierr); 584 ierr = MatDestroy(&B);CHKERRQ(ierr); 585 } 586 PetscFunctionReturn(0); 587 } 588 589 590 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 591 { 592 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 593 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 594 PetscInt *ai = a->i,*ailen = a->ilen; 595 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 596 MatScalar *ap,*aa = a->a; 597 598 PetscFunctionBegin; 599 for (k=0; k<m; k++) { /* loop over rows */ 600 row = im[k]; brow = row/bs; 601 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 602 if (row >= A->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 603 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 604 nrow = ailen[brow]; 605 for (l=0; l<n; l++) { /* loop over columns */ 606 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 607 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 608 col = in[l]; 609 bcol = col/bs; 610 cidx = col%bs; 611 ridx = row%bs; 612 high = nrow; 613 low = 0; /* assume unsorted */ 614 while (high-low > 5) { 615 t = (low+high)/2; 616 if (rp[t] > bcol) high = t; 617 else low = t; 618 } 619 for (i=low; i<high; i++) { 620 if (rp[i] > bcol) break; 621 if (rp[i] == bcol) { 622 *v++ = ap[bs2*i+bs*cidx+ridx]; 623 goto finished; 624 } 625 } 626 *v++ = 0.0; 627 finished:; 628 } 629 } 630 PetscFunctionReturn(0); 631 } 632 633 PetscErrorCode MatPermute_SeqSBAIJ(Mat A,IS rowp,IS colp,Mat *B) 634 { 635 Mat C; 636 PetscErrorCode ierr; 637 638 PetscFunctionBegin; 639 ierr = MatConvert(A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 640 ierr = MatPermute(C,rowp,colp,B);CHKERRQ(ierr); 641 ierr = MatDestroy(&C);CHKERRQ(ierr); 642 if (rowp == colp) { 643 ierr = MatConvert(*B,MATSEQSBAIJ,MAT_INPLACE_MATRIX,B);CHKERRQ(ierr); 644 } 645 PetscFunctionReturn(0); 646 } 647 648 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 649 { 650 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 651 PetscErrorCode ierr; 652 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 653 PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen; 654 PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 655 PetscBool roworiented=a->roworiented; 656 const PetscScalar *value = v; 657 MatScalar *ap,*aa = a->a,*bap; 658 659 PetscFunctionBegin; 660 if (roworiented) stepval = (n-1)*bs; 661 else stepval = (m-1)*bs; 662 663 for (k=0; k<m; k++) { /* loop over added rows */ 664 row = im[k]; 665 if (row < 0) continue; 666 if (PetscUnlikelyDebug(row >= a->mbs)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index row too large %D max %D",row,a->mbs-1); 667 rp = aj + ai[row]; 668 ap = aa + bs2*ai[row]; 669 rmax = imax[row]; 670 nrow = ailen[row]; 671 low = 0; 672 high = nrow; 673 for (l=0; l<n; l++) { /* loop over added columns */ 674 if (in[l] < 0) continue; 675 col = in[l]; 676 if (PetscUnlikelyDebug(col >= a->nbs)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Block index column too large %D max %D",col,a->nbs-1); 677 if (col < row) { 678 if (a->ignore_ltriangular) continue; /* ignore lower triangular block */ 679 else SETERRQ(PETSC_COMM_SELF,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,PETSC_TRUE)"); 680 } 681 if (roworiented) value = v + k*(stepval+bs)*bs + l*bs; 682 else value = v + l*(stepval+bs)*bs + k*bs; 683 684 if (col <= lastcol) low = 0; 685 else high = nrow; 686 687 lastcol = col; 688 while (high-low > 7) { 689 t = (low+high)/2; 690 if (rp[t] > col) high = t; 691 else low = t; 692 } 693 for (i=low; i<high; i++) { 694 if (rp[i] > col) break; 695 if (rp[i] == col) { 696 bap = ap + bs2*i; 697 if (roworiented) { 698 if (is == ADD_VALUES) { 699 for (ii=0; ii<bs; ii++,value+=stepval) { 700 for (jj=ii; jj<bs2; jj+=bs) { 701 bap[jj] += *value++; 702 } 703 } 704 } else { 705 for (ii=0; ii<bs; ii++,value+=stepval) { 706 for (jj=ii; jj<bs2; jj+=bs) { 707 bap[jj] = *value++; 708 } 709 } 710 } 711 } else { 712 if (is == ADD_VALUES) { 713 for (ii=0; ii<bs; ii++,value+=stepval) { 714 for (jj=0; jj<bs; jj++) { 715 *bap++ += *value++; 716 } 717 } 718 } else { 719 for (ii=0; ii<bs; ii++,value+=stepval) { 720 for (jj=0; jj<bs; jj++) { 721 *bap++ = *value++; 722 } 723 } 724 } 725 } 726 goto noinsert2; 727 } 728 } 729 if (nonew == 1) goto noinsert2; 730 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new block index nonzero block (%D, %D) in the matrix", row, col); 731 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 732 N = nrow++ - 1; high++; 733 /* shift up all the later entries in this row */ 734 ierr = PetscArraymove(rp+i+1,rp+i,N-i+1);CHKERRQ(ierr); 735 ierr = PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));CHKERRQ(ierr); 736 ierr = PetscArrayzero(ap+bs2*i,bs2);CHKERRQ(ierr); 737 rp[i] = col; 738 bap = ap + bs2*i; 739 if (roworiented) { 740 for (ii=0; ii<bs; ii++,value+=stepval) { 741 for (jj=ii; jj<bs2; jj+=bs) { 742 bap[jj] = *value++; 743 } 744 } 745 } else { 746 for (ii=0; ii<bs; ii++,value+=stepval) { 747 for (jj=0; jj<bs; jj++) { 748 *bap++ = *value++; 749 } 750 } 751 } 752 noinsert2:; 753 low = i; 754 } 755 ailen[row] = nrow; 756 } 757 PetscFunctionReturn(0); 758 } 759 760 /* 761 This is not yet used 762 */ 763 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 764 { 765 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 766 PetscErrorCode ierr; 767 const PetscInt *ai = a->i, *aj = a->j,*cols; 768 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 769 PetscBool flag; 770 771 PetscFunctionBegin; 772 ierr = PetscMalloc1(m,&ns);CHKERRQ(ierr); 773 while (i < m) { 774 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 775 /* Limits the number of elements in a node to 'a->inode.limit' */ 776 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 777 nzy = ai[j+1] - ai[j]; 778 if (nzy != (nzx - j + i)) break; 779 ierr = PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag);CHKERRQ(ierr); 780 if (!flag) break; 781 } 782 ns[node_count++] = blk_size; 783 784 i = j; 785 } 786 if (!a->inode.size && m && node_count > .9*m) { 787 ierr = PetscFree(ns);CHKERRQ(ierr); 788 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 789 } else { 790 a->inode.node_count = node_count; 791 792 ierr = PetscMalloc1(node_count,&a->inode.size);CHKERRQ(ierr); 793 ierr = PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt));CHKERRQ(ierr); 794 ierr = PetscArraycpy(a->inode.size,ns,node_count);CHKERRQ(ierr); 795 ierr = PetscFree(ns);CHKERRQ(ierr); 796 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 797 798 /* count collections of adjacent columns in each inode */ 799 row = 0; 800 cnt = 0; 801 for (i=0; i<node_count; i++) { 802 cols = aj + ai[row] + a->inode.size[i]; 803 nz = ai[row+1] - ai[row] - a->inode.size[i]; 804 for (j=1; j<nz; j++) { 805 if (cols[j] != cols[j-1]+1) cnt++; 806 } 807 cnt++; 808 row += a->inode.size[i]; 809 } 810 ierr = PetscMalloc1(2*cnt,&counts);CHKERRQ(ierr); 811 cnt = 0; 812 row = 0; 813 for (i=0; i<node_count; i++) { 814 cols = aj + ai[row] + a->inode.size[i]; 815 counts[2*cnt] = cols[0]; 816 nz = ai[row+1] - ai[row] - a->inode.size[i]; 817 cnt2 = 1; 818 for (j=1; j<nz; j++) { 819 if (cols[j] != cols[j-1]+1) { 820 counts[2*(cnt++)+1] = cnt2; 821 counts[2*cnt] = cols[j]; 822 cnt2 = 1; 823 } else cnt2++; 824 } 825 counts[2*(cnt++)+1] = cnt2; 826 row += a->inode.size[i]; 827 } 828 ierr = PetscIntView(2*cnt,counts,NULL);CHKERRQ(ierr); 829 } 830 PetscFunctionReturn(0); 831 } 832 833 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 834 { 835 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 836 PetscErrorCode ierr; 837 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 838 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 839 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 840 MatScalar *aa = a->a,*ap; 841 842 PetscFunctionBegin; 843 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 844 845 if (m) rmax = ailen[0]; 846 for (i=1; i<mbs; i++) { 847 /* move each row back by the amount of empty slots (fshift) before it*/ 848 fshift += imax[i-1] - ailen[i-1]; 849 rmax = PetscMax(rmax,ailen[i]); 850 if (fshift) { 851 ip = aj + ai[i]; 852 ap = aa + bs2*ai[i]; 853 N = ailen[i]; 854 ierr = PetscArraymove(ip-fshift,ip,N);CHKERRQ(ierr); 855 ierr = PetscArraymove(ap-bs2*fshift,ap,bs2*N);CHKERRQ(ierr); 856 } 857 ai[i] = ai[i-1] + ailen[i-1]; 858 } 859 if (mbs) { 860 fshift += imax[mbs-1] - ailen[mbs-1]; 861 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 862 } 863 /* reset ilen and imax for each row */ 864 for (i=0; i<mbs; i++) { 865 ailen[i] = imax[i] = ai[i+1] - ai[i]; 866 } 867 a->nz = ai[mbs]; 868 869 /* diagonals may have moved, reset it */ 870 if (a->diag) { 871 ierr = PetscArraycpy(a->diag,ai,mbs);CHKERRQ(ierr); 872 } 873 if (fshift && a->nounused == -1) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D block size %D, %D unneeded", m, A->cmap->n, A->rmap->bs, fshift*bs2); 874 875 ierr = PetscInfo5(A,"Matrix size: %D X %D, block size %D; storage space: %D unneeded, %D used\n",m,A->rmap->N,A->rmap->bs,fshift*bs2,a->nz*bs2);CHKERRQ(ierr); 876 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 877 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 878 879 A->info.mallocs += a->reallocs; 880 a->reallocs = 0; 881 A->info.nz_unneeded = (PetscReal)fshift*bs2; 882 a->idiagvalid = PETSC_FALSE; 883 a->rmax = rmax; 884 885 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 886 if (a->jshort && a->free_jshort) { 887 /* when matrix data structure is changed, previous jshort must be replaced */ 888 ierr = PetscFree(a->jshort);CHKERRQ(ierr); 889 } 890 ierr = PetscMalloc1(a->i[A->rmap->n],&a->jshort);CHKERRQ(ierr); 891 ierr = PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short));CHKERRQ(ierr); 892 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 893 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 894 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 895 a->free_jshort = PETSC_TRUE; 896 } 897 PetscFunctionReturn(0); 898 } 899 900 /* 901 This function returns an array of flags which indicate the locations of contiguous 902 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 903 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 904 Assume: sizes should be long enough to hold all the values. 905 */ 906 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 907 { 908 PetscInt i,j,k,row; 909 PetscBool flg; 910 911 PetscFunctionBegin; 912 for (i=0,j=0; i<n; j++) { 913 row = idx[i]; 914 if (row%bs!=0) { /* Not the begining of a block */ 915 sizes[j] = 1; 916 i++; 917 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 918 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 919 i++; 920 } else { /* Begining of the block, so check if the complete block exists */ 921 flg = PETSC_TRUE; 922 for (k=1; k<bs; k++) { 923 if (row+k != idx[i+k]) { /* break in the block */ 924 flg = PETSC_FALSE; 925 break; 926 } 927 } 928 if (flg) { /* No break in the bs */ 929 sizes[j] = bs; 930 i += bs; 931 } else { 932 sizes[j] = 1; 933 i++; 934 } 935 } 936 } 937 *bs_max = j; 938 PetscFunctionReturn(0); 939 } 940 941 942 /* Only add/insert a(i,j) with i<=j (blocks). 943 Any a(i,j) with i>j input by user is ingored. 944 */ 945 946 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 947 { 948 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 949 PetscErrorCode ierr; 950 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 951 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 952 PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 953 PetscInt ridx,cidx,bs2=a->bs2; 954 MatScalar *ap,value,*aa=a->a,*bap; 955 956 PetscFunctionBegin; 957 for (k=0; k<m; k++) { /* loop over added rows */ 958 row = im[k]; /* row number */ 959 brow = row/bs; /* block row number */ 960 if (row < 0) continue; 961 if (PetscUnlikelyDebug(row >= A->rmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 962 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 963 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 964 rmax = imax[brow]; /* maximum space allocated for this row */ 965 nrow = ailen[brow]; /* actual length of this row */ 966 low = 0; 967 high = nrow; 968 for (l=0; l<n; l++) { /* loop over added columns */ 969 if (in[l] < 0) continue; 970 if (PetscUnlikelyDebug(in[l] >= A->cmap->N)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->N-1); 971 col = in[l]; 972 bcol = col/bs; /* block col number */ 973 974 if (brow > bcol) { 975 if (a->ignore_ltriangular) continue; /* ignore lower triangular values */ 976 else SETERRQ(PETSC_COMM_SELF,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,PETSC_TRUE)"); 977 } 978 979 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 980 if ((brow==bcol && ridx<=cidx) || (brow<bcol)) { 981 /* element value a(k,l) */ 982 if (roworiented) value = v[l + k*n]; 983 else value = v[k + l*m]; 984 985 /* move pointer bap to a(k,l) quickly and add/insert value */ 986 if (col <= lastcol) low = 0; 987 else high = nrow; 988 989 lastcol = col; 990 while (high-low > 7) { 991 t = (low+high)/2; 992 if (rp[t] > bcol) high = t; 993 else low = t; 994 } 995 for (i=low; i<high; i++) { 996 if (rp[i] > bcol) break; 997 if (rp[i] == bcol) { 998 bap = ap + bs2*i + bs*cidx + ridx; 999 if (is == ADD_VALUES) *bap += value; 1000 else *bap = value; 1001 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1002 if (brow == bcol && ridx < cidx) { 1003 bap = ap + bs2*i + bs*ridx + cidx; 1004 if (is == ADD_VALUES) *bap += value; 1005 else *bap = value; 1006 } 1007 goto noinsert1; 1008 } 1009 } 1010 1011 if (nonew == 1) goto noinsert1; 1012 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 1013 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 1014 1015 N = nrow++ - 1; high++; 1016 /* shift up all the later entries in this row */ 1017 ierr = PetscArraymove(rp+i+1,rp+i,N-i+1);CHKERRQ(ierr); 1018 ierr = PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));CHKERRQ(ierr); 1019 ierr = PetscArrayzero(ap+bs2*i,bs2);CHKERRQ(ierr); 1020 rp[i] = bcol; 1021 ap[bs2*i + bs*cidx + ridx] = value; 1022 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1023 if (brow == bcol && ridx < cidx) { 1024 ap[bs2*i + bs*ridx + cidx] = value; 1025 } 1026 A->nonzerostate++; 1027 noinsert1:; 1028 low = i; 1029 } 1030 } /* end of loop over added columns */ 1031 ailen[brow] = nrow; 1032 } /* end of loop over added rows */ 1033 PetscFunctionReturn(0); 1034 } 1035 1036 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1037 { 1038 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1039 Mat outA; 1040 PetscErrorCode ierr; 1041 PetscBool row_identity; 1042 1043 PetscFunctionBegin; 1044 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1045 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1046 if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1047 if (inA->rmap->bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix block size %D is not supported",inA->rmap->bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */ 1048 1049 outA = inA; 1050 inA->factortype = MAT_FACTOR_ICC; 1051 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 1052 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 1053 1054 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1055 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr); 1056 1057 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1058 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1059 a->row = row; 1060 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1061 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1062 a->col = row; 1063 1064 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1065 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1066 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 1067 1068 if (!a->solve_work) { 1069 ierr = PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);CHKERRQ(ierr); 1070 ierr = PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1071 } 1072 1073 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1074 PetscFunctionReturn(0); 1075 } 1076 1077 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1078 { 1079 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data; 1080 PetscInt i,nz,n; 1081 PetscErrorCode ierr; 1082 1083 PetscFunctionBegin; 1084 nz = baij->maxnz; 1085 n = mat->cmap->n; 1086 for (i=0; i<nz; i++) baij->j[i] = indices[i]; 1087 1088 baij->nz = nz; 1089 for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i]; 1090 1091 ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1092 PetscFunctionReturn(0); 1093 } 1094 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(), and the columns indices MUST be sorted. 1112 1113 MUST be called before any calls to MatSetValues() 1114 1115 .seealso: MatCreateSeqSBAIJ 1116 @*/ 1117 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1118 { 1119 PetscErrorCode ierr; 1120 1121 PetscFunctionBegin; 1122 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1123 PetscValidPointer(indices,2); 1124 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 1125 PetscFunctionReturn(0); 1126 } 1127 1128 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1129 { 1130 PetscErrorCode ierr; 1131 PetscBool isbaij; 1132 1133 PetscFunctionBegin; 1134 ierr = PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"");CHKERRQ(ierr); 1135 if (!isbaij) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name); 1136 /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */ 1137 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 1138 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1139 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1140 1141 if (a->i[a->mbs] != b->i[b->mbs]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1142 if (a->mbs != b->mbs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different"); 1143 if (a->bs2 != b->bs2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size"); 1144 ierr = PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs]);CHKERRQ(ierr); 1145 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 1146 } else { 1147 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1148 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1149 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1150 } 1151 PetscFunctionReturn(0); 1152 } 1153 1154 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A) 1155 { 1156 PetscErrorCode ierr; 1157 1158 PetscFunctionBegin; 1159 ierr = MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL);CHKERRQ(ierr); 1160 PetscFunctionReturn(0); 1161 } 1162 1163 static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1164 { 1165 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1166 1167 PetscFunctionBegin; 1168 *array = a->a; 1169 PetscFunctionReturn(0); 1170 } 1171 1172 static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1173 { 1174 PetscFunctionBegin; 1175 *array = NULL; 1176 PetscFunctionReturn(0); 1177 } 1178 1179 PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz) 1180 { 1181 PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs; 1182 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data; 1183 Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data; 1184 PetscErrorCode ierr; 1185 1186 PetscFunctionBegin; 1187 /* Set the number of nonzeros in the new matrix */ 1188 ierr = MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 1189 PetscFunctionReturn(0); 1190 } 1191 1192 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1193 { 1194 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data; 1195 PetscErrorCode ierr; 1196 PetscInt bs=Y->rmap->bs,bs2=bs*bs; 1197 PetscBLASInt one = 1; 1198 1199 PetscFunctionBegin; 1200 if (str == SAME_NONZERO_PATTERN) { 1201 PetscScalar alpha = a; 1202 PetscBLASInt bnz; 1203 ierr = PetscBLASIntCast(x->nz*bs2,&bnz);CHKERRQ(ierr); 1204 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 1205 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 1206 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1207 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);CHKERRQ(ierr); 1208 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1209 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);CHKERRQ(ierr); 1210 } else { 1211 Mat B; 1212 PetscInt *nnz; 1213 if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size"); 1214 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1215 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 1216 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 1217 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 1218 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 1219 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 1220 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 1221 ierr = MatSetType(B,((PetscObject)Y)->type_name);CHKERRQ(ierr); 1222 ierr = MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);CHKERRQ(ierr); 1223 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1224 1225 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 1226 1227 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 1228 ierr = PetscFree(nnz);CHKERRQ(ierr); 1229 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1230 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 1231 } 1232 PetscFunctionReturn(0); 1233 } 1234 1235 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1236 { 1237 PetscFunctionBegin; 1238 *flg = PETSC_TRUE; 1239 PetscFunctionReturn(0); 1240 } 1241 1242 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1243 { 1244 PetscFunctionBegin; 1245 *flg = PETSC_TRUE; 1246 PetscFunctionReturn(0); 1247 } 1248 1249 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1250 { 1251 PetscFunctionBegin; 1252 *flg = PETSC_FALSE; 1253 PetscFunctionReturn(0); 1254 } 1255 1256 PetscErrorCode MatConjugate_SeqSBAIJ(Mat A) 1257 { 1258 #if defined(PETSC_USE_COMPLEX) 1259 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1260 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1261 MatScalar *aa = a->a; 1262 1263 PetscFunctionBegin; 1264 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1265 #else 1266 PetscFunctionBegin; 1267 #endif 1268 PetscFunctionReturn(0); 1269 } 1270 1271 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1272 { 1273 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1274 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1275 MatScalar *aa = a->a; 1276 1277 PetscFunctionBegin; 1278 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1279 PetscFunctionReturn(0); 1280 } 1281 1282 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1283 { 1284 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1285 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1286 MatScalar *aa = a->a; 1287 1288 PetscFunctionBegin; 1289 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1290 PetscFunctionReturn(0); 1291 } 1292 1293 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1294 { 1295 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1296 PetscErrorCode ierr; 1297 PetscInt i,j,k,count; 1298 PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col; 1299 PetscScalar zero = 0.0; 1300 MatScalar *aa; 1301 const PetscScalar *xx; 1302 PetscScalar *bb; 1303 PetscBool *zeroed,vecs = PETSC_FALSE; 1304 1305 PetscFunctionBegin; 1306 /* fix right hand side if needed */ 1307 if (x && b) { 1308 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1309 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1310 vecs = PETSC_TRUE; 1311 } 1312 1313 /* zero the columns */ 1314 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1315 for (i=0; i<is_n; i++) { 1316 if (is_idx[i] < 0 || is_idx[i] >= A->rmap->N) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range",is_idx[i]); 1317 zeroed[is_idx[i]] = PETSC_TRUE; 1318 } 1319 if (vecs) { 1320 for (i=0; i<A->rmap->N; i++) { 1321 row = i/bs; 1322 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1323 for (k=0; k<bs; k++) { 1324 col = bs*baij->j[j] + k; 1325 if (col <= i) continue; 1326 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1327 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col]; 1328 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i]; 1329 } 1330 } 1331 } 1332 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]]; 1333 } 1334 1335 for (i=0; i<A->rmap->N; i++) { 1336 if (!zeroed[i]) { 1337 row = i/bs; 1338 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1339 for (k=0; k<bs; k++) { 1340 col = bs*baij->j[j] + k; 1341 if (zeroed[col]) { 1342 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1343 aa[0] = 0.0; 1344 } 1345 } 1346 } 1347 } 1348 } 1349 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1350 if (vecs) { 1351 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1352 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1353 } 1354 1355 /* zero the rows */ 1356 for (i=0; i<is_n; i++) { 1357 row = is_idx[i]; 1358 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1359 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1360 for (k=0; k<count; k++) { 1361 aa[0] = zero; 1362 aa += bs; 1363 } 1364 if (diag != 0.0) { 1365 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1366 } 1367 } 1368 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1369 PetscFunctionReturn(0); 1370 } 1371 1372 PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a) 1373 { 1374 PetscErrorCode ierr; 1375 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data; 1376 1377 PetscFunctionBegin; 1378 if (!Y->preallocated || !aij->nz) { 1379 ierr = MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);CHKERRQ(ierr); 1380 } 1381 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 1382 PetscFunctionReturn(0); 1383 } 1384 1385 /* -------------------------------------------------------------------*/ 1386 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1387 MatGetRow_SeqSBAIJ, 1388 MatRestoreRow_SeqSBAIJ, 1389 MatMult_SeqSBAIJ_N, 1390 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1391 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1392 MatMultAdd_SeqSBAIJ_N, 1393 NULL, 1394 NULL, 1395 NULL, 1396 /* 10*/ NULL, 1397 NULL, 1398 MatCholeskyFactor_SeqSBAIJ, 1399 MatSOR_SeqSBAIJ, 1400 MatTranspose_SeqSBAIJ, 1401 /* 15*/ MatGetInfo_SeqSBAIJ, 1402 MatEqual_SeqSBAIJ, 1403 MatGetDiagonal_SeqSBAIJ, 1404 MatDiagonalScale_SeqSBAIJ, 1405 MatNorm_SeqSBAIJ, 1406 /* 20*/ NULL, 1407 MatAssemblyEnd_SeqSBAIJ, 1408 MatSetOption_SeqSBAIJ, 1409 MatZeroEntries_SeqSBAIJ, 1410 /* 24*/ NULL, 1411 NULL, 1412 NULL, 1413 NULL, 1414 NULL, 1415 /* 29*/ MatSetUp_SeqSBAIJ, 1416 NULL, 1417 NULL, 1418 NULL, 1419 NULL, 1420 /* 34*/ MatDuplicate_SeqSBAIJ, 1421 NULL, 1422 NULL, 1423 NULL, 1424 MatICCFactor_SeqSBAIJ, 1425 /* 39*/ MatAXPY_SeqSBAIJ, 1426 MatCreateSubMatrices_SeqSBAIJ, 1427 MatIncreaseOverlap_SeqSBAIJ, 1428 MatGetValues_SeqSBAIJ, 1429 MatCopy_SeqSBAIJ, 1430 /* 44*/ NULL, 1431 MatScale_SeqSBAIJ, 1432 MatShift_SeqSBAIJ, 1433 NULL, 1434 MatZeroRowsColumns_SeqSBAIJ, 1435 /* 49*/ NULL, 1436 MatGetRowIJ_SeqSBAIJ, 1437 MatRestoreRowIJ_SeqSBAIJ, 1438 NULL, 1439 NULL, 1440 /* 54*/ NULL, 1441 NULL, 1442 NULL, 1443 MatPermute_SeqSBAIJ, 1444 MatSetValuesBlocked_SeqSBAIJ, 1445 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1446 NULL, 1447 NULL, 1448 NULL, 1449 NULL, 1450 /* 64*/ NULL, 1451 NULL, 1452 NULL, 1453 NULL, 1454 NULL, 1455 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1456 NULL, 1457 MatConvert_MPISBAIJ_Basic, 1458 NULL, 1459 NULL, 1460 /* 74*/ NULL, 1461 NULL, 1462 NULL, 1463 NULL, 1464 NULL, 1465 /* 79*/ NULL, 1466 NULL, 1467 NULL, 1468 MatGetInertia_SeqSBAIJ, 1469 MatLoad_SeqSBAIJ, 1470 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1471 MatIsHermitian_SeqSBAIJ, 1472 MatIsStructurallySymmetric_SeqSBAIJ, 1473 NULL, 1474 NULL, 1475 /* 89*/ NULL, 1476 NULL, 1477 NULL, 1478 NULL, 1479 NULL, 1480 /* 94*/ NULL, 1481 NULL, 1482 NULL, 1483 NULL, 1484 NULL, 1485 /* 99*/ NULL, 1486 NULL, 1487 NULL, 1488 MatConjugate_SeqSBAIJ, 1489 NULL, 1490 /*104*/ NULL, 1491 MatRealPart_SeqSBAIJ, 1492 MatImaginaryPart_SeqSBAIJ, 1493 MatGetRowUpperTriangular_SeqSBAIJ, 1494 MatRestoreRowUpperTriangular_SeqSBAIJ, 1495 /*109*/ NULL, 1496 NULL, 1497 NULL, 1498 NULL, 1499 MatMissingDiagonal_SeqSBAIJ, 1500 /*114*/ NULL, 1501 NULL, 1502 NULL, 1503 NULL, 1504 NULL, 1505 /*119*/ NULL, 1506 NULL, 1507 NULL, 1508 NULL, 1509 NULL, 1510 /*124*/ NULL, 1511 NULL, 1512 NULL, 1513 NULL, 1514 NULL, 1515 /*129*/ NULL, 1516 NULL, 1517 NULL, 1518 NULL, 1519 NULL, 1520 /*134*/ NULL, 1521 NULL, 1522 NULL, 1523 NULL, 1524 NULL, 1525 /*139*/ MatSetBlockSizes_Default, 1526 NULL, 1527 NULL, 1528 NULL, 1529 NULL, 1530 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ 1531 }; 1532 1533 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1534 { 1535 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1536 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1537 PetscErrorCode ierr; 1538 1539 PetscFunctionBegin; 1540 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1541 1542 /* allocate space for values if not already there */ 1543 if (!aij->saved_values) { 1544 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 1545 } 1546 1547 /* copy values over */ 1548 ierr = PetscArraycpy(aij->saved_values,aij->a,nz);CHKERRQ(ierr); 1549 PetscFunctionReturn(0); 1550 } 1551 1552 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1553 { 1554 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1555 PetscErrorCode ierr; 1556 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1557 1558 PetscFunctionBegin; 1559 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1560 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1561 1562 /* copy values over */ 1563 ierr = PetscArraycpy(aij->a,aij->saved_values,nz);CHKERRQ(ierr); 1564 PetscFunctionReturn(0); 1565 } 1566 1567 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1568 { 1569 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1570 PetscErrorCode ierr; 1571 PetscInt i,mbs,nbs,bs2; 1572 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1573 1574 PetscFunctionBegin; 1575 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1576 1577 ierr = MatSetBlockSize(B,PetscAbs(bs));CHKERRQ(ierr); 1578 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1579 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1580 if (B->rmap->N > B->cmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"SEQSBAIJ matrix cannot have more rows %D than columns %D",B->rmap->N,B->cmap->N); 1581 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1582 1583 B->preallocated = PETSC_TRUE; 1584 1585 mbs = B->rmap->N/bs; 1586 nbs = B->cmap->n/bs; 1587 bs2 = bs*bs; 1588 1589 if (mbs*bs != B->rmap->N || nbs*bs!=B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1590 1591 if (nz == MAT_SKIP_ALLOCATION) { 1592 skipallocation = PETSC_TRUE; 1593 nz = 0; 1594 } 1595 1596 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1597 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1598 if (nnz) { 1599 for (i=0; i<mbs; i++) { 1600 if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1601 if (nnz[i] > nbs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %D value %D block rowlength %D",i,nnz[i],nbs); 1602 } 1603 } 1604 1605 B->ops->mult = MatMult_SeqSBAIJ_N; 1606 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1607 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1608 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1609 1610 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1611 if (!flg) { 1612 switch (bs) { 1613 case 1: 1614 B->ops->mult = MatMult_SeqSBAIJ_1; 1615 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1616 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1617 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1618 break; 1619 case 2: 1620 B->ops->mult = MatMult_SeqSBAIJ_2; 1621 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1622 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1623 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1624 break; 1625 case 3: 1626 B->ops->mult = MatMult_SeqSBAIJ_3; 1627 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1628 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1629 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1630 break; 1631 case 4: 1632 B->ops->mult = MatMult_SeqSBAIJ_4; 1633 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1634 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1635 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1636 break; 1637 case 5: 1638 B->ops->mult = MatMult_SeqSBAIJ_5; 1639 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1640 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1641 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1642 break; 1643 case 6: 1644 B->ops->mult = MatMult_SeqSBAIJ_6; 1645 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1646 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1647 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1648 break; 1649 case 7: 1650 B->ops->mult = MatMult_SeqSBAIJ_7; 1651 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1652 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1653 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1654 break; 1655 } 1656 } 1657 1658 b->mbs = mbs; 1659 b->nbs = nbs; 1660 if (!skipallocation) { 1661 if (!b->imax) { 1662 ierr = PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);CHKERRQ(ierr); 1663 1664 b->free_imax_ilen = PETSC_TRUE; 1665 1666 ierr = PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1667 } 1668 if (!nnz) { 1669 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1670 else if (nz <= 0) nz = 1; 1671 nz = PetscMin(nbs,nz); 1672 for (i=0; i<mbs; i++) b->imax[i] = nz; 1673 ierr = PetscIntMultError(nz,mbs,&nz);CHKERRQ(ierr); 1674 } else { 1675 PetscInt64 nz64 = 0; 1676 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 1677 ierr = PetscIntCast(nz64,&nz);CHKERRQ(ierr); 1678 } 1679 /* b->ilen will count nonzeros in each block row so far. */ 1680 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1681 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1682 1683 /* allocate the matrix space */ 1684 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1685 ierr = PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);CHKERRQ(ierr); 1686 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1687 ierr = PetscArrayzero(b->a,nz*bs2);CHKERRQ(ierr); 1688 ierr = PetscArrayzero(b->j,nz);CHKERRQ(ierr); 1689 1690 b->singlemalloc = PETSC_TRUE; 1691 1692 /* pointer to beginning of each row */ 1693 b->i[0] = 0; 1694 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1695 1696 b->free_a = PETSC_TRUE; 1697 b->free_ij = PETSC_TRUE; 1698 } else { 1699 b->free_a = PETSC_FALSE; 1700 b->free_ij = PETSC_FALSE; 1701 } 1702 1703 b->bs2 = bs2; 1704 b->nz = 0; 1705 b->maxnz = nz; 1706 b->inew = NULL; 1707 b->jnew = NULL; 1708 b->anew = NULL; 1709 b->a2anew = NULL; 1710 b->permute = PETSC_FALSE; 1711 1712 B->was_assembled = PETSC_FALSE; 1713 B->assembled = PETSC_FALSE; 1714 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1715 PetscFunctionReturn(0); 1716 } 1717 1718 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1719 { 1720 PetscInt i,j,m,nz,anz, nz_max=0,*nnz; 1721 PetscScalar *values=NULL; 1722 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1723 PetscErrorCode ierr; 1724 1725 PetscFunctionBegin; 1726 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs); 1727 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1728 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1729 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1730 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1731 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1732 m = B->rmap->n/bs; 1733 1734 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]); 1735 ierr = PetscMalloc1(m+1,&nnz);CHKERRQ(ierr); 1736 for (i=0; i<m; i++) { 1737 nz = ii[i+1] - ii[i]; 1738 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz); 1739 anz = 0; 1740 for (j=0; j<nz; j++) { 1741 /* count only values on the diagonal or above */ 1742 if (jj[ii[i] + j] >= i) { 1743 anz = nz - j; 1744 break; 1745 } 1746 } 1747 nz_max = PetscMax(nz_max,anz); 1748 nnz[i] = anz; 1749 } 1750 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1751 ierr = PetscFree(nnz);CHKERRQ(ierr); 1752 1753 values = (PetscScalar*)V; 1754 if (!values) { 1755 ierr = PetscCalloc1(bs*bs*nz_max,&values);CHKERRQ(ierr); 1756 } 1757 for (i=0; i<m; i++) { 1758 PetscInt ncols = ii[i+1] - ii[i]; 1759 const PetscInt *icols = jj + ii[i]; 1760 if (!roworiented || bs == 1) { 1761 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1762 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr); 1763 } else { 1764 for (j=0; j<ncols; j++) { 1765 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1766 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr); 1767 } 1768 } 1769 } 1770 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); } 1771 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1772 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1773 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1774 PetscFunctionReturn(0); 1775 } 1776 1777 /* 1778 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1779 */ 1780 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1781 { 1782 PetscErrorCode ierr; 1783 PetscBool flg = PETSC_FALSE; 1784 PetscInt bs = B->rmap->bs; 1785 1786 PetscFunctionBegin; 1787 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1788 if (flg) bs = 8; 1789 1790 if (!natural) { 1791 switch (bs) { 1792 case 1: 1793 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1794 break; 1795 case 2: 1796 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1797 break; 1798 case 3: 1799 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1800 break; 1801 case 4: 1802 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1803 break; 1804 case 5: 1805 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1806 break; 1807 case 6: 1808 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1809 break; 1810 case 7: 1811 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1812 break; 1813 default: 1814 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1815 break; 1816 } 1817 } else { 1818 switch (bs) { 1819 case 1: 1820 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1821 break; 1822 case 2: 1823 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1824 break; 1825 case 3: 1826 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1827 break; 1828 case 4: 1829 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1830 break; 1831 case 5: 1832 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1833 break; 1834 case 6: 1835 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1836 break; 1837 case 7: 1838 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1839 break; 1840 default: 1841 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1842 break; 1843 } 1844 } 1845 PetscFunctionReturn(0); 1846 } 1847 1848 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*); 1849 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat,MatType,MatReuse,Mat*); 1850 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type) 1851 { 1852 PetscFunctionBegin; 1853 *type = MATSOLVERPETSC; 1854 PetscFunctionReturn(0); 1855 } 1856 1857 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1858 { 1859 PetscInt n = A->rmap->n; 1860 PetscErrorCode ierr; 1861 1862 PetscFunctionBegin; 1863 #if defined(PETSC_USE_COMPLEX) 1864 if (A->hermitian && !A->symmetric && (ftype == MAT_FACTOR_CHOLESKY||ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian CHOLESKY or ICC Factor is not supported"); 1865 #endif 1866 1867 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 1868 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1869 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1870 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1871 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1872 1873 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1874 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1875 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);CHKERRQ(ierr); 1876 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);CHKERRQ(ierr); 1877 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1878 1879 (*B)->factortype = ftype; 1880 (*B)->canuseordering = PETSC_TRUE; 1881 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 1882 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 1883 ierr = PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc);CHKERRQ(ierr); 1884 PetscFunctionReturn(0); 1885 } 1886 1887 /*@C 1888 MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored 1889 1890 Not Collective 1891 1892 Input Parameter: 1893 . mat - a MATSEQSBAIJ matrix 1894 1895 Output Parameter: 1896 . array - pointer to the data 1897 1898 Level: intermediate 1899 1900 .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray() 1901 @*/ 1902 PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array) 1903 { 1904 PetscErrorCode ierr; 1905 1906 PetscFunctionBegin; 1907 ierr = PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1908 PetscFunctionReturn(0); 1909 } 1910 1911 /*@C 1912 MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray() 1913 1914 Not Collective 1915 1916 Input Parameters: 1917 + mat - a MATSEQSBAIJ matrix 1918 - array - pointer to the data 1919 1920 Level: intermediate 1921 1922 .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray() 1923 @*/ 1924 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array) 1925 { 1926 PetscErrorCode ierr; 1927 1928 PetscFunctionBegin; 1929 ierr = PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1930 PetscFunctionReturn(0); 1931 } 1932 1933 /*MC 1934 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1935 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1936 1937 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1938 can call MatSetOption(Mat, MAT_HERMITIAN). 1939 1940 Options Database Keys: 1941 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1942 1943 Notes: 1944 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1945 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1946 the options database -mat_ignore_lower_triangular false it will generate an error if you try to set a value in the lower triangular portion. 1947 1948 The number of rows in the matrix must be less than or equal to the number of columns 1949 1950 Level: beginner 1951 1952 .seealso: MatCreateSeqSBAIJ(), MatType, MATMPISBAIJ 1953 M*/ 1954 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1955 { 1956 Mat_SeqSBAIJ *b; 1957 PetscErrorCode ierr; 1958 PetscMPIInt size; 1959 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1960 1961 PetscFunctionBegin; 1962 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRMPI(ierr); 1963 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1964 1965 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1966 B->data = (void*)b; 1967 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1968 1969 B->ops->destroy = MatDestroy_SeqSBAIJ; 1970 B->ops->view = MatView_SeqSBAIJ; 1971 b->row = NULL; 1972 b->icol = NULL; 1973 b->reallocs = 0; 1974 b->saved_values = NULL; 1975 b->inode.limit = 5; 1976 b->inode.max_limit = 5; 1977 1978 b->roworiented = PETSC_TRUE; 1979 b->nonew = 0; 1980 b->diag = NULL; 1981 b->solve_work = NULL; 1982 b->mult_work = NULL; 1983 B->spptr = NULL; 1984 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1985 b->keepnonzeropattern = PETSC_FALSE; 1986 1987 b->inew = NULL; 1988 b->jnew = NULL; 1989 b->anew = NULL; 1990 b->a2anew = NULL; 1991 b->permute = PETSC_FALSE; 1992 1993 b->ignore_ltriangular = PETSC_TRUE; 1994 1995 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1996 1997 b->getrow_utriangular = PETSC_FALSE; 1998 1999 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 2000 2001 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);CHKERRQ(ierr); 2002 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);CHKERRQ(ierr); 2003 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 2004 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 2005 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 2006 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 2007 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 2008 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 2009 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr); 2010 #if defined(PETSC_HAVE_ELEMENTAL) 2011 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);CHKERRQ(ierr); 2012 #endif 2013 #if defined(PETSC_HAVE_SCALAPACK) 2014 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK);CHKERRQ(ierr); 2015 #endif 2016 2017 B->symmetric = PETSC_TRUE; 2018 B->structurally_symmetric = PETSC_TRUE; 2019 B->symmetric_set = PETSC_TRUE; 2020 B->structurally_symmetric_set = PETSC_TRUE; 2021 B->symmetric_eternal = PETSC_TRUE; 2022 #if defined(PETSC_USE_COMPLEX) 2023 B->hermitian = PETSC_FALSE; 2024 B->hermitian_set = PETSC_FALSE; 2025 #else 2026 B->hermitian = PETSC_TRUE; 2027 B->hermitian_set = PETSC_TRUE; 2028 #endif 2029 2030 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 2031 2032 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 2033 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 2034 if (no_unroll) { 2035 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 2036 } 2037 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 2038 if (no_inode) { 2039 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 2040 } 2041 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 2042 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2043 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 2044 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 2045 PetscFunctionReturn(0); 2046 } 2047 2048 /*@C 2049 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2050 compressed row) format. For good matrix assembly performance the 2051 user should preallocate the matrix storage by setting the parameter nz 2052 (or the array nnz). By setting these parameters accurately, performance 2053 during matrix assembly can be increased by more than a factor of 50. 2054 2055 Collective on Mat 2056 2057 Input Parameters: 2058 + B - the symmetric matrix 2059 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2060 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2061 . nz - number of block nonzeros per block row (same for all rows) 2062 - nnz - array containing the number of block nonzeros in the upper triangular plus 2063 diagonal portion of each block (possibly different for each block row) or NULL 2064 2065 Options Database Keys: 2066 + -mat_no_unroll - uses code that does not unroll the loops in the 2067 block calculations (much slower) 2068 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2069 2070 Level: intermediate 2071 2072 Notes: 2073 Specify the preallocated storage with either nz or nnz (not both). 2074 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2075 allocation. See Users-Manual: ch_mat for details. 2076 2077 You can call MatGetInfo() to get information on how effective the preallocation was; 2078 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2079 You can also run with the option -info and look for messages with the string 2080 malloc in them to see if additional memory allocation was needed. 2081 2082 If the nnz parameter is given then the nz parameter is ignored 2083 2084 2085 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2086 @*/ 2087 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2088 { 2089 PetscErrorCode ierr; 2090 2091 PetscFunctionBegin; 2092 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2093 PetscValidType(B,1); 2094 PetscValidLogicalCollectiveInt(B,bs,2); 2095 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2096 PetscFunctionReturn(0); 2097 } 2098 2099 /*@C 2100 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values 2101 2102 Input Parameters: 2103 + B - the matrix 2104 . bs - size of block, the blocks are ALWAYS square. 2105 . i - the indices into j for the start of each local row (starts with zero) 2106 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2107 - v - optional values in the matrix 2108 2109 Level: advanced 2110 2111 Notes: 2112 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2113 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2114 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2115 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2116 block column and the second index is over columns within a block. 2117 2118 Any entries below the diagonal are ignored 2119 2120 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries 2121 and usually the numerical values as well 2122 2123 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ 2124 @*/ 2125 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2126 { 2127 PetscErrorCode ierr; 2128 2129 PetscFunctionBegin; 2130 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2131 PetscValidType(B,1); 2132 PetscValidLogicalCollectiveInt(B,bs,2); 2133 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2134 PetscFunctionReturn(0); 2135 } 2136 2137 /*@C 2138 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2139 compressed row) format. For good matrix assembly performance the 2140 user should preallocate the matrix storage by setting the parameter nz 2141 (or the array nnz). By setting these parameters accurately, performance 2142 during matrix assembly can be increased by more than a factor of 50. 2143 2144 Collective 2145 2146 Input Parameters: 2147 + comm - MPI communicator, set to PETSC_COMM_SELF 2148 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2149 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2150 . m - number of rows, or number of columns 2151 . nz - number of block nonzeros per block row (same for all rows) 2152 - nnz - array containing the number of block nonzeros in the upper triangular plus 2153 diagonal portion of each block (possibly different for each block row) or NULL 2154 2155 Output Parameter: 2156 . A - the symmetric matrix 2157 2158 Options Database Keys: 2159 + -mat_no_unroll - uses code that does not unroll the loops in the 2160 block calculations (much slower) 2161 - -mat_block_size - size of the blocks to use 2162 2163 Level: intermediate 2164 2165 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2166 MatXXXXSetPreallocation() paradigm instead of this routine directly. 2167 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2168 2169 Notes: 2170 The number of rows and columns must be divisible by blocksize. 2171 This matrix type does not support complex Hermitian operation. 2172 2173 Specify the preallocated storage with either nz or nnz (not both). 2174 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2175 allocation. See Users-Manual: ch_mat for details. 2176 2177 If the nnz parameter is given then the nz parameter is ignored 2178 2179 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2180 @*/ 2181 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2182 { 2183 PetscErrorCode ierr; 2184 2185 PetscFunctionBegin; 2186 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2187 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2188 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2189 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2190 PetscFunctionReturn(0); 2191 } 2192 2193 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2194 { 2195 Mat C; 2196 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2197 PetscErrorCode ierr; 2198 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2199 2200 PetscFunctionBegin; 2201 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2202 2203 *B = NULL; 2204 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2205 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2206 ierr = MatSetBlockSizesFromMats(C,A,A);CHKERRQ(ierr); 2207 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2208 c = (Mat_SeqSBAIJ*)C->data; 2209 2210 C->preallocated = PETSC_TRUE; 2211 C->factortype = A->factortype; 2212 c->row = NULL; 2213 c->icol = NULL; 2214 c->saved_values = NULL; 2215 c->keepnonzeropattern = a->keepnonzeropattern; 2216 C->assembled = PETSC_TRUE; 2217 2218 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2219 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2220 c->bs2 = a->bs2; 2221 c->mbs = a->mbs; 2222 c->nbs = a->nbs; 2223 2224 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2225 c->imax = a->imax; 2226 c->ilen = a->ilen; 2227 c->free_imax_ilen = PETSC_FALSE; 2228 } else { 2229 ierr = PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);CHKERRQ(ierr); 2230 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2231 for (i=0; i<mbs; i++) { 2232 c->imax[i] = a->imax[i]; 2233 c->ilen[i] = a->ilen[i]; 2234 } 2235 c->free_imax_ilen = PETSC_TRUE; 2236 } 2237 2238 /* allocate the matrix space */ 2239 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2240 ierr = PetscMalloc1(bs2*nz,&c->a);CHKERRQ(ierr); 2241 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2242 c->i = a->i; 2243 c->j = a->j; 2244 c->singlemalloc = PETSC_FALSE; 2245 c->free_a = PETSC_TRUE; 2246 c->free_ij = PETSC_FALSE; 2247 c->parent = A; 2248 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2249 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2250 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2251 } else { 2252 ierr = PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);CHKERRQ(ierr); 2253 ierr = PetscArraycpy(c->i,a->i,mbs+1);CHKERRQ(ierr); 2254 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2255 c->singlemalloc = PETSC_TRUE; 2256 c->free_a = PETSC_TRUE; 2257 c->free_ij = PETSC_TRUE; 2258 } 2259 if (mbs > 0) { 2260 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2261 ierr = PetscArraycpy(c->j,a->j,nz);CHKERRQ(ierr); 2262 } 2263 if (cpvalues == MAT_COPY_VALUES) { 2264 ierr = PetscArraycpy(c->a,a->a,bs2*nz);CHKERRQ(ierr); 2265 } else { 2266 ierr = PetscArrayzero(c->a,bs2*nz);CHKERRQ(ierr); 2267 } 2268 if (a->jshort) { 2269 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2270 /* if the parent matrix is reassembled, this child matrix will never notice */ 2271 ierr = PetscMalloc1(nz,&c->jshort);CHKERRQ(ierr); 2272 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2273 ierr = PetscArraycpy(c->jshort,a->jshort,nz);CHKERRQ(ierr); 2274 2275 c->free_jshort = PETSC_TRUE; 2276 } 2277 } 2278 2279 c->roworiented = a->roworiented; 2280 c->nonew = a->nonew; 2281 2282 if (a->diag) { 2283 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2284 c->diag = a->diag; 2285 c->free_diag = PETSC_FALSE; 2286 } else { 2287 ierr = PetscMalloc1(mbs,&c->diag);CHKERRQ(ierr); 2288 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2289 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2290 c->free_diag = PETSC_TRUE; 2291 } 2292 } 2293 c->nz = a->nz; 2294 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2295 c->solve_work = NULL; 2296 c->mult_work = NULL; 2297 2298 *B = C; 2299 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2300 PetscFunctionReturn(0); 2301 } 2302 2303 /* Used for both SeqBAIJ and SeqSBAIJ matrices */ 2304 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary 2305 2306 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer) 2307 { 2308 PetscErrorCode ierr; 2309 PetscBool isbinary; 2310 2311 PetscFunctionBegin; 2312 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 2313 if (!isbinary) SETERRQ2(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)mat)->type_name); 2314 ierr = MatLoad_SeqSBAIJ_Binary(mat,viewer);CHKERRQ(ierr); 2315 PetscFunctionReturn(0); 2316 } 2317 2318 /*@ 2319 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2320 (upper triangular entries in CSR format) provided by the user. 2321 2322 Collective 2323 2324 Input Parameters: 2325 + comm - must be an MPI communicator of size 1 2326 . bs - size of block 2327 . m - number of rows 2328 . n - number of columns 2329 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of block elements in that row block row of the matrix 2330 . j - column indices 2331 - a - matrix values 2332 2333 Output Parameter: 2334 . mat - the matrix 2335 2336 Level: advanced 2337 2338 Notes: 2339 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2340 once the matrix is destroyed 2341 2342 You cannot set new nonzero locations into this matrix, that will generate an error. 2343 2344 The i and j indices are 0 based 2345 2346 When block size is greater than 1 the matrix values must be stored using the SBAIJ storage format (see the SBAIJ code to determine this). For block size of 1 2347 it is the regular CSR format excluding the lower triangular elements. 2348 2349 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2350 2351 @*/ 2352 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 2353 { 2354 PetscErrorCode ierr; 2355 PetscInt ii; 2356 Mat_SeqSBAIJ *sbaij; 2357 2358 PetscFunctionBegin; 2359 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2360 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2361 2362 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2363 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2364 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2365 ierr = MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 2366 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2367 ierr = PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);CHKERRQ(ierr); 2368 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2369 2370 sbaij->i = i; 2371 sbaij->j = j; 2372 sbaij->a = a; 2373 2374 sbaij->singlemalloc = PETSC_FALSE; 2375 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2376 sbaij->free_a = PETSC_FALSE; 2377 sbaij->free_ij = PETSC_FALSE; 2378 sbaij->free_imax_ilen = PETSC_TRUE; 2379 2380 for (ii=0; ii<m; ii++) { 2381 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2382 if (PetscUnlikelyDebug(i[ii+1] - i[ii] < 0)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2383 } 2384 if (PetscDefined(USE_DEBUG)) { 2385 for (ii=0; ii<sbaij->i[m]; ii++) { 2386 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2387 if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2388 } 2389 } 2390 2391 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2392 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2393 PetscFunctionReturn(0); 2394 } 2395 2396 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2397 { 2398 PetscErrorCode ierr; 2399 PetscMPIInt size; 2400 2401 PetscFunctionBegin; 2402 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 2403 if (size == 1 && scall == MAT_REUSE_MATRIX) { 2404 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2405 } else { 2406 ierr = MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2407 } 2408 PetscFunctionReturn(0); 2409 } 2410