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 (nz) *nz = 0; 313 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);} 314 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);} 315 PetscFunctionReturn(0); 316 } 317 318 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 319 { 320 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 321 322 PetscFunctionBegin; 323 a->getrow_utriangular = PETSC_TRUE; 324 PetscFunctionReturn(0); 325 } 326 327 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) 328 { 329 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 330 331 PetscFunctionBegin; 332 a->getrow_utriangular = PETSC_FALSE; 333 PetscFunctionReturn(0); 334 } 335 336 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B) 337 { 338 PetscErrorCode ierr; 339 340 PetscFunctionBegin; 341 if (reuse == MAT_INITIAL_MATRIX) { 342 ierr = MatDuplicate(A,MAT_COPY_VALUES,B);CHKERRQ(ierr); 343 } else if (reuse == MAT_REUSE_MATRIX) { 344 ierr = MatCopy(A,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 345 } 346 PetscFunctionReturn(0); 347 } 348 349 PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 350 { 351 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 352 PetscErrorCode ierr; 353 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2; 354 PetscViewerFormat format; 355 PetscInt *diag; 356 357 PetscFunctionBegin; 358 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 359 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 360 ierr = PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);CHKERRQ(ierr); 361 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 362 Mat aij; 363 const char *matname; 364 365 if (A->factortype && bs>1) { 366 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); 367 PetscFunctionReturn(0); 368 } 369 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);CHKERRQ(ierr); 370 ierr = PetscObjectGetName((PetscObject)A,&matname);CHKERRQ(ierr); 371 ierr = PetscObjectSetName((PetscObject)aij,matname);CHKERRQ(ierr); 372 ierr = MatView(aij,viewer);CHKERRQ(ierr); 373 ierr = MatDestroy(&aij);CHKERRQ(ierr); 374 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 375 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 376 for (i=0; i<a->mbs; i++) { 377 for (j=0; j<bs; j++) { 378 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 379 for (k=a->i[i]; k<a->i[i+1]; k++) { 380 for (l=0; l<bs; l++) { 381 #if defined(PETSC_USE_COMPLEX) 382 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 383 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l, 384 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 385 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 386 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l, 387 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 388 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 389 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 390 } 391 #else 392 if (a->a[bs2*k + l*bs + j] != 0.0) { 393 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 394 } 395 #endif 396 } 397 } 398 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 399 } 400 } 401 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 402 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 403 PetscFunctionReturn(0); 404 } else { 405 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 406 if (A->factortype) { /* for factored matrix */ 407 if (bs>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet"); 408 409 diag=a->diag; 410 for (i=0; i<a->mbs; i++) { /* for row block i */ 411 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 412 /* diagonal entry */ 413 #if defined(PETSC_USE_COMPLEX) 414 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) { 415 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); 416 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) { 417 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); 418 } else { 419 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 420 } 421 #else 422 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]));CHKERRQ(ierr); 423 #endif 424 /* off-diagonal entries */ 425 for (k=a->i[i]; k<a->i[i+1]-1; k++) { 426 #if defined(PETSC_USE_COMPLEX) 427 if (PetscImaginaryPart(a->a[k]) > 0.0) { 428 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 429 } else if (PetscImaginaryPart(a->a[k]) < 0.0) { 430 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 431 } else { 432 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]));CHKERRQ(ierr); 433 } 434 #else 435 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[k],(double)a->a[k]);CHKERRQ(ierr); 436 #endif 437 } 438 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 439 } 440 441 } else { /* for non-factored matrix */ 442 for (i=0; i<a->mbs; i++) { /* for row block i */ 443 for (j=0; j<bs; j++) { /* for row bs*i + j */ 444 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 445 for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */ 446 for (l=0; l<bs; l++) { /* for column */ 447 #if defined(PETSC_USE_COMPLEX) 448 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 449 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l, 450 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 451 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 452 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l, 453 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 454 } else { 455 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 456 } 457 #else 458 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 459 #endif 460 } 461 } 462 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 463 } 464 } 465 } 466 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 467 } 468 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 469 PetscFunctionReturn(0); 470 } 471 472 #include <petscdraw.h> 473 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 474 { 475 Mat A = (Mat) Aa; 476 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 477 PetscErrorCode ierr; 478 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 479 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 480 MatScalar *aa; 481 PetscViewer viewer; 482 483 PetscFunctionBegin; 484 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 485 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 486 487 /* loop over matrix elements drawing boxes */ 488 489 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 490 ierr = PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric");CHKERRQ(ierr); 491 /* Blue for negative, Cyan for zero and Red for positive */ 492 color = PETSC_DRAW_BLUE; 493 for (i=0,row=0; i<mbs; i++,row+=bs) { 494 for (j=a->i[i]; j<a->i[i+1]; j++) { 495 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 496 x_l = a->j[j]*bs; x_r = x_l + 1.0; 497 aa = a->a + j*bs2; 498 for (k=0; k<bs; k++) { 499 for (l=0; l<bs; l++) { 500 if (PetscRealPart(*aa++) >= 0.) continue; 501 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 502 } 503 } 504 } 505 } 506 color = PETSC_DRAW_CYAN; 507 for (i=0,row=0; i<mbs; i++,row+=bs) { 508 for (j=a->i[i]; j<a->i[i+1]; j++) { 509 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 510 x_l = a->j[j]*bs; x_r = x_l + 1.0; 511 aa = a->a + j*bs2; 512 for (k=0; k<bs; k++) { 513 for (l=0; l<bs; l++) { 514 if (PetscRealPart(*aa++) != 0.) continue; 515 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 516 } 517 } 518 } 519 } 520 color = PETSC_DRAW_RED; 521 for (i=0,row=0; i<mbs; i++,row+=bs) { 522 for (j=a->i[i]; j<a->i[i+1]; j++) { 523 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 524 x_l = a->j[j]*bs; x_r = x_l + 1.0; 525 aa = a->a + j*bs2; 526 for (k=0; k<bs; k++) { 527 for (l=0; l<bs; l++) { 528 if (PetscRealPart(*aa++) <= 0.) continue; 529 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 530 } 531 } 532 } 533 } 534 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 535 PetscFunctionReturn(0); 536 } 537 538 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 539 { 540 PetscErrorCode ierr; 541 PetscReal xl,yl,xr,yr,w,h; 542 PetscDraw draw; 543 PetscBool isnull; 544 545 PetscFunctionBegin; 546 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 547 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 548 if (isnull) PetscFunctionReturn(0); 549 550 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 551 xr += w; yr += h; xl = -w; yl = -h; 552 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 553 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 554 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 555 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 556 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 557 PetscFunctionReturn(0); 558 } 559 560 /* Used for both MPIBAIJ and MPISBAIJ matrices */ 561 #define MatView_SeqSBAIJ_Binary MatView_SeqBAIJ_Binary 562 563 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 564 { 565 PetscErrorCode ierr; 566 PetscBool iascii,isbinary,isdraw; 567 568 PetscFunctionBegin; 569 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 570 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 571 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 572 if (iascii) { 573 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 574 } else if (isbinary) { 575 ierr = MatView_SeqSBAIJ_Binary(A,viewer);CHKERRQ(ierr); 576 } else if (isdraw) { 577 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 578 } else { 579 Mat B; 580 const char *matname; 581 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 582 ierr = PetscObjectGetName((PetscObject)A,&matname);CHKERRQ(ierr); 583 ierr = PetscObjectSetName((PetscObject)B,matname);CHKERRQ(ierr); 584 ierr = MatView(B,viewer);CHKERRQ(ierr); 585 ierr = MatDestroy(&B);CHKERRQ(ierr); 586 } 587 PetscFunctionReturn(0); 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] corresponding 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 beginning 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 /* Only add/insert a(i,j) with i<=j (blocks). 942 Any a(i,j) with i>j input by user is ingored. 943 */ 944 945 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 946 { 947 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 948 PetscErrorCode ierr; 949 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 950 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 951 PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 952 PetscInt ridx,cidx,bs2=a->bs2; 953 MatScalar *ap,value,*aa=a->a,*bap; 954 955 PetscFunctionBegin; 956 for (k=0; k<m; k++) { /* loop over added rows */ 957 row = im[k]; /* row number */ 958 brow = row/bs; /* block row number */ 959 if (row < 0) continue; 960 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); 961 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 962 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 963 rmax = imax[brow]; /* maximum space allocated for this row */ 964 nrow = ailen[brow]; /* actual length of this row */ 965 low = 0; 966 high = nrow; 967 for (l=0; l<n; l++) { /* loop over added columns */ 968 if (in[l] < 0) continue; 969 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); 970 col = in[l]; 971 bcol = col/bs; /* block col number */ 972 973 if (brow > bcol) { 974 if (a->ignore_ltriangular) continue; /* ignore lower triangular values */ 975 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)"); 976 } 977 978 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 979 if ((brow==bcol && ridx<=cidx) || (brow<bcol)) { 980 /* element value a(k,l) */ 981 if (roworiented) value = v[l + k*n]; 982 else value = v[k + l*m]; 983 984 /* move pointer bap to a(k,l) quickly and add/insert value */ 985 if (col <= lastcol) low = 0; 986 else high = nrow; 987 988 lastcol = col; 989 while (high-low > 7) { 990 t = (low+high)/2; 991 if (rp[t] > bcol) high = t; 992 else low = t; 993 } 994 for (i=low; i<high; i++) { 995 if (rp[i] > bcol) break; 996 if (rp[i] == bcol) { 997 bap = ap + bs2*i + bs*cidx + ridx; 998 if (is == ADD_VALUES) *bap += value; 999 else *bap = value; 1000 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1001 if (brow == bcol && ridx < cidx) { 1002 bap = ap + bs2*i + bs*ridx + cidx; 1003 if (is == ADD_VALUES) *bap += value; 1004 else *bap = value; 1005 } 1006 goto noinsert1; 1007 } 1008 } 1009 1010 if (nonew == 1) goto noinsert1; 1011 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 1012 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 1013 1014 N = nrow++ - 1; high++; 1015 /* shift up all the later entries in this row */ 1016 ierr = PetscArraymove(rp+i+1,rp+i,N-i+1);CHKERRQ(ierr); 1017 ierr = PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));CHKERRQ(ierr); 1018 ierr = PetscArrayzero(ap+bs2*i,bs2);CHKERRQ(ierr); 1019 rp[i] = bcol; 1020 ap[bs2*i + bs*cidx + ridx] = value; 1021 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1022 if (brow == bcol && ridx < cidx) { 1023 ap[bs2*i + bs*ridx + cidx] = value; 1024 } 1025 A->nonzerostate++; 1026 noinsert1:; 1027 low = i; 1028 } 1029 } /* end of loop over added columns */ 1030 ailen[brow] = nrow; 1031 } /* end of loop over added rows */ 1032 PetscFunctionReturn(0); 1033 } 1034 1035 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1036 { 1037 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1038 Mat outA; 1039 PetscErrorCode ierr; 1040 PetscBool row_identity; 1041 1042 PetscFunctionBegin; 1043 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1044 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1045 if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1046 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()! */ 1047 1048 outA = inA; 1049 inA->factortype = MAT_FACTOR_ICC; 1050 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 1051 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 1052 1053 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1054 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr); 1055 1056 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1057 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1058 a->row = row; 1059 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1060 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1061 a->col = row; 1062 1063 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1064 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1065 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 1066 1067 if (!a->solve_work) { 1068 ierr = PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);CHKERRQ(ierr); 1069 ierr = PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1070 } 1071 1072 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1073 PetscFunctionReturn(0); 1074 } 1075 1076 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1077 { 1078 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data; 1079 PetscInt i,nz,n; 1080 PetscErrorCode ierr; 1081 1082 PetscFunctionBegin; 1083 nz = baij->maxnz; 1084 n = mat->cmap->n; 1085 for (i=0; i<nz; i++) baij->j[i] = indices[i]; 1086 1087 baij->nz = nz; 1088 for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i]; 1089 1090 ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1091 PetscFunctionReturn(0); 1092 } 1093 1094 /*@ 1095 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1096 in the matrix. 1097 1098 Input Parameters: 1099 + mat - the SeqSBAIJ matrix 1100 - indices - the column indices 1101 1102 Level: advanced 1103 1104 Notes: 1105 This can be called if you have precomputed the nonzero structure of the 1106 matrix and want to provide it to the matrix object to improve the performance 1107 of the MatSetValues() operation. 1108 1109 You MUST have set the correct numbers of nonzeros per row in the call to 1110 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1111 1112 MUST be called before any calls to MatSetValues() 1113 1114 .seealso: MatCreateSeqSBAIJ 1115 @*/ 1116 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1117 { 1118 PetscErrorCode ierr; 1119 1120 PetscFunctionBegin; 1121 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1122 PetscValidPointer(indices,2); 1123 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 1124 PetscFunctionReturn(0); 1125 } 1126 1127 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1128 { 1129 PetscErrorCode ierr; 1130 PetscBool isbaij; 1131 1132 PetscFunctionBegin; 1133 ierr = PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"");CHKERRQ(ierr); 1134 if (!isbaij) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name); 1135 /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */ 1136 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 1137 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1138 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1139 1140 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"); 1141 if (a->mbs != b->mbs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different"); 1142 if (a->bs2 != b->bs2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size"); 1143 ierr = PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs]);CHKERRQ(ierr); 1144 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 1145 } else { 1146 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1147 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1148 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1149 } 1150 PetscFunctionReturn(0); 1151 } 1152 1153 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A) 1154 { 1155 PetscErrorCode ierr; 1156 1157 PetscFunctionBegin; 1158 ierr = MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL);CHKERRQ(ierr); 1159 PetscFunctionReturn(0); 1160 } 1161 1162 static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1163 { 1164 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1165 1166 PetscFunctionBegin; 1167 *array = a->a; 1168 PetscFunctionReturn(0); 1169 } 1170 1171 static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1172 { 1173 PetscFunctionBegin; 1174 *array = NULL; 1175 PetscFunctionReturn(0); 1176 } 1177 1178 PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz) 1179 { 1180 PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs; 1181 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data; 1182 Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data; 1183 PetscErrorCode ierr; 1184 1185 PetscFunctionBegin; 1186 /* Set the number of nonzeros in the new matrix */ 1187 ierr = MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 1188 PetscFunctionReturn(0); 1189 } 1190 1191 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1192 { 1193 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data; 1194 PetscErrorCode ierr; 1195 PetscInt bs=Y->rmap->bs,bs2=bs*bs; 1196 PetscBLASInt one = 1; 1197 1198 PetscFunctionBegin; 1199 if (str == SAME_NONZERO_PATTERN) { 1200 PetscScalar alpha = a; 1201 PetscBLASInt bnz; 1202 ierr = PetscBLASIntCast(x->nz*bs2,&bnz);CHKERRQ(ierr); 1203 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 1204 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 1205 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1206 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);CHKERRQ(ierr); 1207 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1208 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);CHKERRQ(ierr); 1209 } else { 1210 Mat B; 1211 PetscInt *nnz; 1212 if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size"); 1213 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1214 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 1215 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 1216 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 1217 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 1218 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 1219 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 1220 ierr = MatSetType(B,((PetscObject)Y)->type_name);CHKERRQ(ierr); 1221 ierr = MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);CHKERRQ(ierr); 1222 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1223 1224 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 1225 1226 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 1227 ierr = PetscFree(nnz);CHKERRQ(ierr); 1228 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1229 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 1230 } 1231 PetscFunctionReturn(0); 1232 } 1233 1234 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1235 { 1236 PetscFunctionBegin; 1237 *flg = PETSC_TRUE; 1238 PetscFunctionReturn(0); 1239 } 1240 1241 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1242 { 1243 PetscFunctionBegin; 1244 *flg = PETSC_TRUE; 1245 PetscFunctionReturn(0); 1246 } 1247 1248 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1249 { 1250 PetscFunctionBegin; 1251 *flg = PETSC_FALSE; 1252 PetscFunctionReturn(0); 1253 } 1254 1255 PetscErrorCode MatConjugate_SeqSBAIJ(Mat A) 1256 { 1257 #if defined(PETSC_USE_COMPLEX) 1258 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1259 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1260 MatScalar *aa = a->a; 1261 1262 PetscFunctionBegin; 1263 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1264 #else 1265 PetscFunctionBegin; 1266 #endif 1267 PetscFunctionReturn(0); 1268 } 1269 1270 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1271 { 1272 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1273 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1274 MatScalar *aa = a->a; 1275 1276 PetscFunctionBegin; 1277 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1278 PetscFunctionReturn(0); 1279 } 1280 1281 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1282 { 1283 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1284 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1285 MatScalar *aa = a->a; 1286 1287 PetscFunctionBegin; 1288 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1289 PetscFunctionReturn(0); 1290 } 1291 1292 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1293 { 1294 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1295 PetscErrorCode ierr; 1296 PetscInt i,j,k,count; 1297 PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col; 1298 PetscScalar zero = 0.0; 1299 MatScalar *aa; 1300 const PetscScalar *xx; 1301 PetscScalar *bb; 1302 PetscBool *zeroed,vecs = PETSC_FALSE; 1303 1304 PetscFunctionBegin; 1305 /* fix right hand side if needed */ 1306 if (x && b) { 1307 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1308 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1309 vecs = PETSC_TRUE; 1310 } 1311 1312 /* zero the columns */ 1313 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1314 for (i=0; i<is_n; i++) { 1315 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]); 1316 zeroed[is_idx[i]] = PETSC_TRUE; 1317 } 1318 if (vecs) { 1319 for (i=0; i<A->rmap->N; i++) { 1320 row = i/bs; 1321 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1322 for (k=0; k<bs; k++) { 1323 col = bs*baij->j[j] + k; 1324 if (col <= i) continue; 1325 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1326 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col]; 1327 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i]; 1328 } 1329 } 1330 } 1331 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]]; 1332 } 1333 1334 for (i=0; i<A->rmap->N; i++) { 1335 if (!zeroed[i]) { 1336 row = i/bs; 1337 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1338 for (k=0; k<bs; k++) { 1339 col = bs*baij->j[j] + k; 1340 if (zeroed[col]) { 1341 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1342 aa[0] = 0.0; 1343 } 1344 } 1345 } 1346 } 1347 } 1348 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1349 if (vecs) { 1350 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1351 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1352 } 1353 1354 /* zero the rows */ 1355 for (i=0; i<is_n; i++) { 1356 row = is_idx[i]; 1357 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1358 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1359 for (k=0; k<count; k++) { 1360 aa[0] = zero; 1361 aa += bs; 1362 } 1363 if (diag != 0.0) { 1364 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1365 } 1366 } 1367 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1368 PetscFunctionReturn(0); 1369 } 1370 1371 PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a) 1372 { 1373 PetscErrorCode ierr; 1374 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data; 1375 1376 PetscFunctionBegin; 1377 if (!Y->preallocated || !aij->nz) { 1378 ierr = MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);CHKERRQ(ierr); 1379 } 1380 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 1381 PetscFunctionReturn(0); 1382 } 1383 1384 /* -------------------------------------------------------------------*/ 1385 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1386 MatGetRow_SeqSBAIJ, 1387 MatRestoreRow_SeqSBAIJ, 1388 MatMult_SeqSBAIJ_N, 1389 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1390 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1391 MatMultAdd_SeqSBAIJ_N, 1392 NULL, 1393 NULL, 1394 NULL, 1395 /* 10*/ NULL, 1396 NULL, 1397 MatCholeskyFactor_SeqSBAIJ, 1398 MatSOR_SeqSBAIJ, 1399 MatTranspose_SeqSBAIJ, 1400 /* 15*/ MatGetInfo_SeqSBAIJ, 1401 MatEqual_SeqSBAIJ, 1402 MatGetDiagonal_SeqSBAIJ, 1403 MatDiagonalScale_SeqSBAIJ, 1404 MatNorm_SeqSBAIJ, 1405 /* 20*/ NULL, 1406 MatAssemblyEnd_SeqSBAIJ, 1407 MatSetOption_SeqSBAIJ, 1408 MatZeroEntries_SeqSBAIJ, 1409 /* 24*/ NULL, 1410 NULL, 1411 NULL, 1412 NULL, 1413 NULL, 1414 /* 29*/ MatSetUp_SeqSBAIJ, 1415 NULL, 1416 NULL, 1417 NULL, 1418 NULL, 1419 /* 34*/ MatDuplicate_SeqSBAIJ, 1420 NULL, 1421 NULL, 1422 NULL, 1423 MatICCFactor_SeqSBAIJ, 1424 /* 39*/ MatAXPY_SeqSBAIJ, 1425 MatCreateSubMatrices_SeqSBAIJ, 1426 MatIncreaseOverlap_SeqSBAIJ, 1427 MatGetValues_SeqSBAIJ, 1428 MatCopy_SeqSBAIJ, 1429 /* 44*/ NULL, 1430 MatScale_SeqSBAIJ, 1431 MatShift_SeqSBAIJ, 1432 NULL, 1433 MatZeroRowsColumns_SeqSBAIJ, 1434 /* 49*/ NULL, 1435 MatGetRowIJ_SeqSBAIJ, 1436 MatRestoreRowIJ_SeqSBAIJ, 1437 NULL, 1438 NULL, 1439 /* 54*/ NULL, 1440 NULL, 1441 NULL, 1442 MatPermute_SeqSBAIJ, 1443 MatSetValuesBlocked_SeqSBAIJ, 1444 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1445 NULL, 1446 NULL, 1447 NULL, 1448 NULL, 1449 /* 64*/ NULL, 1450 NULL, 1451 NULL, 1452 NULL, 1453 NULL, 1454 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1455 NULL, 1456 MatConvert_MPISBAIJ_Basic, 1457 NULL, 1458 NULL, 1459 /* 74*/ NULL, 1460 NULL, 1461 NULL, 1462 NULL, 1463 NULL, 1464 /* 79*/ NULL, 1465 NULL, 1466 NULL, 1467 MatGetInertia_SeqSBAIJ, 1468 MatLoad_SeqSBAIJ, 1469 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1470 MatIsHermitian_SeqSBAIJ, 1471 MatIsStructurallySymmetric_SeqSBAIJ, 1472 NULL, 1473 NULL, 1474 /* 89*/ NULL, 1475 NULL, 1476 NULL, 1477 NULL, 1478 NULL, 1479 /* 94*/ NULL, 1480 NULL, 1481 NULL, 1482 NULL, 1483 NULL, 1484 /* 99*/ NULL, 1485 NULL, 1486 NULL, 1487 MatConjugate_SeqSBAIJ, 1488 NULL, 1489 /*104*/ NULL, 1490 MatRealPart_SeqSBAIJ, 1491 MatImaginaryPart_SeqSBAIJ, 1492 MatGetRowUpperTriangular_SeqSBAIJ, 1493 MatRestoreRowUpperTriangular_SeqSBAIJ, 1494 /*109*/ NULL, 1495 NULL, 1496 NULL, 1497 NULL, 1498 MatMissingDiagonal_SeqSBAIJ, 1499 /*114*/ NULL, 1500 NULL, 1501 NULL, 1502 NULL, 1503 NULL, 1504 /*119*/ NULL, 1505 NULL, 1506 NULL, 1507 NULL, 1508 NULL, 1509 /*124*/ NULL, 1510 NULL, 1511 NULL, 1512 NULL, 1513 NULL, 1514 /*129*/ NULL, 1515 NULL, 1516 NULL, 1517 NULL, 1518 NULL, 1519 /*134*/ NULL, 1520 NULL, 1521 NULL, 1522 NULL, 1523 NULL, 1524 /*139*/ MatSetBlockSizes_Default, 1525 NULL, 1526 NULL, 1527 NULL, 1528 NULL, 1529 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ 1530 }; 1531 1532 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1533 { 1534 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1535 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1536 PetscErrorCode ierr; 1537 1538 PetscFunctionBegin; 1539 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1540 1541 /* allocate space for values if not already there */ 1542 if (!aij->saved_values) { 1543 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 1544 } 1545 1546 /* copy values over */ 1547 ierr = PetscArraycpy(aij->saved_values,aij->a,nz);CHKERRQ(ierr); 1548 PetscFunctionReturn(0); 1549 } 1550 1551 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1552 { 1553 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1554 PetscErrorCode ierr; 1555 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1556 1557 PetscFunctionBegin; 1558 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1559 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1560 1561 /* copy values over */ 1562 ierr = PetscArraycpy(aij->a,aij->saved_values,nz);CHKERRQ(ierr); 1563 PetscFunctionReturn(0); 1564 } 1565 1566 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1567 { 1568 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1569 PetscErrorCode ierr; 1570 PetscInt i,mbs,nbs,bs2; 1571 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1572 1573 PetscFunctionBegin; 1574 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1575 1576 ierr = MatSetBlockSize(B,PetscAbs(bs));CHKERRQ(ierr); 1577 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1578 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1579 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); 1580 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1581 1582 B->preallocated = PETSC_TRUE; 1583 1584 mbs = B->rmap->N/bs; 1585 nbs = B->cmap->n/bs; 1586 bs2 = bs*bs; 1587 1588 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"); 1589 1590 if (nz == MAT_SKIP_ALLOCATION) { 1591 skipallocation = PETSC_TRUE; 1592 nz = 0; 1593 } 1594 1595 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1596 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1597 if (nnz) { 1598 for (i=0; i<mbs; i++) { 1599 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]); 1600 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); 1601 } 1602 } 1603 1604 B->ops->mult = MatMult_SeqSBAIJ_N; 1605 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1606 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1607 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1608 1609 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1610 if (!flg) { 1611 switch (bs) { 1612 case 1: 1613 B->ops->mult = MatMult_SeqSBAIJ_1; 1614 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1615 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1616 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1617 break; 1618 case 2: 1619 B->ops->mult = MatMult_SeqSBAIJ_2; 1620 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1621 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1622 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1623 break; 1624 case 3: 1625 B->ops->mult = MatMult_SeqSBAIJ_3; 1626 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1627 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1628 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1629 break; 1630 case 4: 1631 B->ops->mult = MatMult_SeqSBAIJ_4; 1632 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1633 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1634 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1635 break; 1636 case 5: 1637 B->ops->mult = MatMult_SeqSBAIJ_5; 1638 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1639 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1640 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1641 break; 1642 case 6: 1643 B->ops->mult = MatMult_SeqSBAIJ_6; 1644 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1645 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1646 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1647 break; 1648 case 7: 1649 B->ops->mult = MatMult_SeqSBAIJ_7; 1650 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1651 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1652 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1653 break; 1654 } 1655 } 1656 1657 b->mbs = mbs; 1658 b->nbs = nbs; 1659 if (!skipallocation) { 1660 if (!b->imax) { 1661 ierr = PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);CHKERRQ(ierr); 1662 1663 b->free_imax_ilen = PETSC_TRUE; 1664 1665 ierr = PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1666 } 1667 if (!nnz) { 1668 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1669 else if (nz <= 0) nz = 1; 1670 nz = PetscMin(nbs,nz); 1671 for (i=0; i<mbs; i++) b->imax[i] = nz; 1672 ierr = PetscIntMultError(nz,mbs,&nz);CHKERRQ(ierr); 1673 } else { 1674 PetscInt64 nz64 = 0; 1675 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 1676 ierr = PetscIntCast(nz64,&nz);CHKERRQ(ierr); 1677 } 1678 /* b->ilen will count nonzeros in each block row so far. */ 1679 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1680 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1681 1682 /* allocate the matrix space */ 1683 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1684 ierr = PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);CHKERRQ(ierr); 1685 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1686 ierr = PetscArrayzero(b->a,nz*bs2);CHKERRQ(ierr); 1687 ierr = PetscArrayzero(b->j,nz);CHKERRQ(ierr); 1688 1689 b->singlemalloc = PETSC_TRUE; 1690 1691 /* pointer to beginning of each row */ 1692 b->i[0] = 0; 1693 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1694 1695 b->free_a = PETSC_TRUE; 1696 b->free_ij = PETSC_TRUE; 1697 } else { 1698 b->free_a = PETSC_FALSE; 1699 b->free_ij = PETSC_FALSE; 1700 } 1701 1702 b->bs2 = bs2; 1703 b->nz = 0; 1704 b->maxnz = nz; 1705 b->inew = NULL; 1706 b->jnew = NULL; 1707 b->anew = NULL; 1708 b->a2anew = NULL; 1709 b->permute = PETSC_FALSE; 1710 1711 B->was_assembled = PETSC_FALSE; 1712 B->assembled = PETSC_FALSE; 1713 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1714 PetscFunctionReturn(0); 1715 } 1716 1717 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1718 { 1719 PetscInt i,j,m,nz,anz, nz_max=0,*nnz; 1720 PetscScalar *values=NULL; 1721 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1722 PetscErrorCode ierr; 1723 1724 PetscFunctionBegin; 1725 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs); 1726 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1727 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1728 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1729 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1730 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1731 m = B->rmap->n/bs; 1732 1733 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]); 1734 ierr = PetscMalloc1(m+1,&nnz);CHKERRQ(ierr); 1735 for (i=0; i<m; i++) { 1736 nz = ii[i+1] - ii[i]; 1737 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz); 1738 anz = 0; 1739 for (j=0; j<nz; j++) { 1740 /* count only values on the diagonal or above */ 1741 if (jj[ii[i] + j] >= i) { 1742 anz = nz - j; 1743 break; 1744 } 1745 } 1746 nz_max = PetscMax(nz_max,anz); 1747 nnz[i] = anz; 1748 } 1749 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1750 ierr = PetscFree(nnz);CHKERRQ(ierr); 1751 1752 values = (PetscScalar*)V; 1753 if (!values) { 1754 ierr = PetscCalloc1(bs*bs*nz_max,&values);CHKERRQ(ierr); 1755 } 1756 for (i=0; i<m; i++) { 1757 PetscInt ncols = ii[i+1] - ii[i]; 1758 const PetscInt *icols = jj + ii[i]; 1759 if (!roworiented || bs == 1) { 1760 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1761 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr); 1762 } else { 1763 for (j=0; j<ncols; j++) { 1764 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1765 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr); 1766 } 1767 } 1768 } 1769 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); } 1770 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1771 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1772 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1773 PetscFunctionReturn(0); 1774 } 1775 1776 /* 1777 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1778 */ 1779 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1780 { 1781 PetscErrorCode ierr; 1782 PetscBool flg = PETSC_FALSE; 1783 PetscInt bs = B->rmap->bs; 1784 1785 PetscFunctionBegin; 1786 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1787 if (flg) bs = 8; 1788 1789 if (!natural) { 1790 switch (bs) { 1791 case 1: 1792 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1793 break; 1794 case 2: 1795 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1796 break; 1797 case 3: 1798 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1799 break; 1800 case 4: 1801 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1802 break; 1803 case 5: 1804 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1805 break; 1806 case 6: 1807 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1808 break; 1809 case 7: 1810 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1811 break; 1812 default: 1813 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1814 break; 1815 } 1816 } else { 1817 switch (bs) { 1818 case 1: 1819 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1820 break; 1821 case 2: 1822 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1823 break; 1824 case 3: 1825 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1826 break; 1827 case 4: 1828 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1829 break; 1830 case 5: 1831 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1832 break; 1833 case 6: 1834 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1835 break; 1836 case 7: 1837 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1838 break; 1839 default: 1840 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1841 break; 1842 } 1843 } 1844 PetscFunctionReturn(0); 1845 } 1846 1847 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*); 1848 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat,MatType,MatReuse,Mat*); 1849 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type) 1850 { 1851 PetscFunctionBegin; 1852 *type = MATSOLVERPETSC; 1853 PetscFunctionReturn(0); 1854 } 1855 1856 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1857 { 1858 PetscInt n = A->rmap->n; 1859 PetscErrorCode ierr; 1860 1861 PetscFunctionBegin; 1862 #if defined(PETSC_USE_COMPLEX) 1863 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"); 1864 #endif 1865 1866 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 1867 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1868 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1869 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1870 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1871 1872 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1873 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1874 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);CHKERRQ(ierr); 1875 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);CHKERRQ(ierr); 1876 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1877 1878 (*B)->factortype = ftype; 1879 (*B)->canuseordering = PETSC_TRUE; 1880 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 1881 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 1882 ierr = PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc);CHKERRQ(ierr); 1883 PetscFunctionReturn(0); 1884 } 1885 1886 /*@C 1887 MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored 1888 1889 Not Collective 1890 1891 Input Parameter: 1892 . mat - a MATSEQSBAIJ matrix 1893 1894 Output Parameter: 1895 . array - pointer to the data 1896 1897 Level: intermediate 1898 1899 .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray() 1900 @*/ 1901 PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array) 1902 { 1903 PetscErrorCode ierr; 1904 1905 PetscFunctionBegin; 1906 ierr = PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1907 PetscFunctionReturn(0); 1908 } 1909 1910 /*@C 1911 MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray() 1912 1913 Not Collective 1914 1915 Input Parameters: 1916 + mat - a MATSEQSBAIJ matrix 1917 - array - pointer to the data 1918 1919 Level: intermediate 1920 1921 .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray() 1922 @*/ 1923 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array) 1924 { 1925 PetscErrorCode ierr; 1926 1927 PetscFunctionBegin; 1928 ierr = PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1929 PetscFunctionReturn(0); 1930 } 1931 1932 /*MC 1933 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1934 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1935 1936 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1937 can call MatSetOption(Mat, MAT_HERMITIAN). 1938 1939 Options Database Keys: 1940 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1941 1942 Notes: 1943 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1944 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1945 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. 1946 1947 The number of rows in the matrix must be less than or equal to the number of columns 1948 1949 Level: beginner 1950 1951 .seealso: MatCreateSeqSBAIJ(), MatType, MATMPISBAIJ 1952 M*/ 1953 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1954 { 1955 Mat_SeqSBAIJ *b; 1956 PetscErrorCode ierr; 1957 PetscMPIInt size; 1958 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1959 1960 PetscFunctionBegin; 1961 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRMPI(ierr); 1962 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1963 1964 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1965 B->data = (void*)b; 1966 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1967 1968 B->ops->destroy = MatDestroy_SeqSBAIJ; 1969 B->ops->view = MatView_SeqSBAIJ; 1970 b->row = NULL; 1971 b->icol = NULL; 1972 b->reallocs = 0; 1973 b->saved_values = NULL; 1974 b->inode.limit = 5; 1975 b->inode.max_limit = 5; 1976 1977 b->roworiented = PETSC_TRUE; 1978 b->nonew = 0; 1979 b->diag = NULL; 1980 b->solve_work = NULL; 1981 b->mult_work = NULL; 1982 B->spptr = NULL; 1983 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1984 b->keepnonzeropattern = PETSC_FALSE; 1985 1986 b->inew = NULL; 1987 b->jnew = NULL; 1988 b->anew = NULL; 1989 b->a2anew = NULL; 1990 b->permute = PETSC_FALSE; 1991 1992 b->ignore_ltriangular = PETSC_TRUE; 1993 1994 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1995 1996 b->getrow_utriangular = PETSC_FALSE; 1997 1998 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 1999 2000 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);CHKERRQ(ierr); 2001 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);CHKERRQ(ierr); 2002 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 2003 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 2004 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 2005 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 2006 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 2007 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 2008 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr); 2009 #if defined(PETSC_HAVE_ELEMENTAL) 2010 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);CHKERRQ(ierr); 2011 #endif 2012 #if defined(PETSC_HAVE_SCALAPACK) 2013 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK);CHKERRQ(ierr); 2014 #endif 2015 2016 B->symmetric = PETSC_TRUE; 2017 B->structurally_symmetric = PETSC_TRUE; 2018 B->symmetric_set = PETSC_TRUE; 2019 B->structurally_symmetric_set = PETSC_TRUE; 2020 B->symmetric_eternal = PETSC_TRUE; 2021 #if defined(PETSC_USE_COMPLEX) 2022 B->hermitian = PETSC_FALSE; 2023 B->hermitian_set = PETSC_FALSE; 2024 #else 2025 B->hermitian = PETSC_TRUE; 2026 B->hermitian_set = PETSC_TRUE; 2027 #endif 2028 2029 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 2030 2031 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 2032 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 2033 if (no_unroll) { 2034 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 2035 } 2036 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 2037 if (no_inode) { 2038 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 2039 } 2040 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 2041 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2042 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 2043 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 2044 PetscFunctionReturn(0); 2045 } 2046 2047 /*@C 2048 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2049 compressed row) format. For good matrix assembly performance the 2050 user should preallocate the matrix storage by setting the parameter nz 2051 (or the array nnz). By setting these parameters accurately, performance 2052 during matrix assembly can be increased by more than a factor of 50. 2053 2054 Collective on Mat 2055 2056 Input Parameters: 2057 + B - the symmetric matrix 2058 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2059 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2060 . nz - number of block nonzeros per block row (same for all rows) 2061 - nnz - array containing the number of block nonzeros in the upper triangular plus 2062 diagonal portion of each block (possibly different for each block row) or NULL 2063 2064 Options Database Keys: 2065 + -mat_no_unroll - uses code that does not unroll the loops in the 2066 block calculations (much slower) 2067 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2068 2069 Level: intermediate 2070 2071 Notes: 2072 Specify the preallocated storage with either nz or nnz (not both). 2073 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2074 allocation. See Users-Manual: ch_mat for details. 2075 2076 You can call MatGetInfo() to get information on how effective the preallocation was; 2077 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2078 You can also run with the option -info and look for messages with the string 2079 malloc in them to see if additional memory allocation was needed. 2080 2081 If the nnz parameter is given then the nz parameter is ignored 2082 2083 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2084 @*/ 2085 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2086 { 2087 PetscErrorCode ierr; 2088 2089 PetscFunctionBegin; 2090 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2091 PetscValidType(B,1); 2092 PetscValidLogicalCollectiveInt(B,bs,2); 2093 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2094 PetscFunctionReturn(0); 2095 } 2096 2097 /*@C 2098 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values 2099 2100 Input Parameters: 2101 + B - the matrix 2102 . bs - size of block, the blocks are ALWAYS square. 2103 . i - the indices into j for the start of each local row (starts with zero) 2104 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2105 - v - optional values in the matrix 2106 2107 Level: advanced 2108 2109 Notes: 2110 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2111 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2112 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2113 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2114 block column and the second index is over columns within a block. 2115 2116 Any entries below the diagonal are ignored 2117 2118 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries 2119 and usually the numerical values as well 2120 2121 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ 2122 @*/ 2123 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2124 { 2125 PetscErrorCode ierr; 2126 2127 PetscFunctionBegin; 2128 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2129 PetscValidType(B,1); 2130 PetscValidLogicalCollectiveInt(B,bs,2); 2131 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2132 PetscFunctionReturn(0); 2133 } 2134 2135 /*@C 2136 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2137 compressed row) format. For good matrix assembly performance the 2138 user should preallocate the matrix storage by setting the parameter nz 2139 (or the array nnz). By setting these parameters accurately, performance 2140 during matrix assembly can be increased by more than a factor of 50. 2141 2142 Collective 2143 2144 Input Parameters: 2145 + comm - MPI communicator, set to PETSC_COMM_SELF 2146 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2147 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2148 . m - number of rows, or number of columns 2149 . nz - number of block nonzeros per block row (same for all rows) 2150 - nnz - array containing the number of block nonzeros in the upper triangular plus 2151 diagonal portion of each block (possibly different for each block row) or NULL 2152 2153 Output Parameter: 2154 . A - the symmetric matrix 2155 2156 Options Database Keys: 2157 + -mat_no_unroll - uses code that does not unroll the loops in the 2158 block calculations (much slower) 2159 - -mat_block_size - size of the blocks to use 2160 2161 Level: intermediate 2162 2163 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2164 MatXXXXSetPreallocation() paradigm instead of this routine directly. 2165 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2166 2167 Notes: 2168 The number of rows and columns must be divisible by blocksize. 2169 This matrix type does not support complex Hermitian operation. 2170 2171 Specify the preallocated storage with either nz or nnz (not both). 2172 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2173 allocation. See Users-Manual: ch_mat for details. 2174 2175 If the nnz parameter is given then the nz parameter is ignored 2176 2177 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2178 @*/ 2179 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2180 { 2181 PetscErrorCode ierr; 2182 2183 PetscFunctionBegin; 2184 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2185 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2186 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2187 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2188 PetscFunctionReturn(0); 2189 } 2190 2191 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2192 { 2193 Mat C; 2194 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2195 PetscErrorCode ierr; 2196 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2197 2198 PetscFunctionBegin; 2199 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2200 2201 *B = NULL; 2202 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2203 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2204 ierr = MatSetBlockSizesFromMats(C,A,A);CHKERRQ(ierr); 2205 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2206 c = (Mat_SeqSBAIJ*)C->data; 2207 2208 C->preallocated = PETSC_TRUE; 2209 C->factortype = A->factortype; 2210 c->row = NULL; 2211 c->icol = NULL; 2212 c->saved_values = NULL; 2213 c->keepnonzeropattern = a->keepnonzeropattern; 2214 C->assembled = PETSC_TRUE; 2215 2216 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2217 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2218 c->bs2 = a->bs2; 2219 c->mbs = a->mbs; 2220 c->nbs = a->nbs; 2221 2222 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2223 c->imax = a->imax; 2224 c->ilen = a->ilen; 2225 c->free_imax_ilen = PETSC_FALSE; 2226 } else { 2227 ierr = PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);CHKERRQ(ierr); 2228 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2229 for (i=0; i<mbs; i++) { 2230 c->imax[i] = a->imax[i]; 2231 c->ilen[i] = a->ilen[i]; 2232 } 2233 c->free_imax_ilen = PETSC_TRUE; 2234 } 2235 2236 /* allocate the matrix space */ 2237 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2238 ierr = PetscMalloc1(bs2*nz,&c->a);CHKERRQ(ierr); 2239 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2240 c->i = a->i; 2241 c->j = a->j; 2242 c->singlemalloc = PETSC_FALSE; 2243 c->free_a = PETSC_TRUE; 2244 c->free_ij = PETSC_FALSE; 2245 c->parent = A; 2246 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2247 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2248 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2249 } else { 2250 ierr = PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);CHKERRQ(ierr); 2251 ierr = PetscArraycpy(c->i,a->i,mbs+1);CHKERRQ(ierr); 2252 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2253 c->singlemalloc = PETSC_TRUE; 2254 c->free_a = PETSC_TRUE; 2255 c->free_ij = PETSC_TRUE; 2256 } 2257 if (mbs > 0) { 2258 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2259 ierr = PetscArraycpy(c->j,a->j,nz);CHKERRQ(ierr); 2260 } 2261 if (cpvalues == MAT_COPY_VALUES) { 2262 ierr = PetscArraycpy(c->a,a->a,bs2*nz);CHKERRQ(ierr); 2263 } else { 2264 ierr = PetscArrayzero(c->a,bs2*nz);CHKERRQ(ierr); 2265 } 2266 if (a->jshort) { 2267 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2268 /* if the parent matrix is reassembled, this child matrix will never notice */ 2269 ierr = PetscMalloc1(nz,&c->jshort);CHKERRQ(ierr); 2270 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2271 ierr = PetscArraycpy(c->jshort,a->jshort,nz);CHKERRQ(ierr); 2272 2273 c->free_jshort = PETSC_TRUE; 2274 } 2275 } 2276 2277 c->roworiented = a->roworiented; 2278 c->nonew = a->nonew; 2279 2280 if (a->diag) { 2281 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2282 c->diag = a->diag; 2283 c->free_diag = PETSC_FALSE; 2284 } else { 2285 ierr = PetscMalloc1(mbs,&c->diag);CHKERRQ(ierr); 2286 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2287 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2288 c->free_diag = PETSC_TRUE; 2289 } 2290 } 2291 c->nz = a->nz; 2292 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2293 c->solve_work = NULL; 2294 c->mult_work = NULL; 2295 2296 *B = C; 2297 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2298 PetscFunctionReturn(0); 2299 } 2300 2301 /* Used for both SeqBAIJ and SeqSBAIJ matrices */ 2302 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary 2303 2304 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer) 2305 { 2306 PetscErrorCode ierr; 2307 PetscBool isbinary; 2308 2309 PetscFunctionBegin; 2310 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 2311 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); 2312 ierr = MatLoad_SeqSBAIJ_Binary(mat,viewer);CHKERRQ(ierr); 2313 PetscFunctionReturn(0); 2314 } 2315 2316 /*@ 2317 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2318 (upper triangular entries in CSR format) provided by the user. 2319 2320 Collective 2321 2322 Input Parameters: 2323 + comm - must be an MPI communicator of size 1 2324 . bs - size of block 2325 . m - number of rows 2326 . n - number of columns 2327 . 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 2328 . j - column indices 2329 - a - matrix values 2330 2331 Output Parameter: 2332 . mat - the matrix 2333 2334 Level: advanced 2335 2336 Notes: 2337 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2338 once the matrix is destroyed 2339 2340 You cannot set new nonzero locations into this matrix, that will generate an error. 2341 2342 The i and j indices are 0 based 2343 2344 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 2345 it is the regular CSR format excluding the lower triangular elements. 2346 2347 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2348 2349 @*/ 2350 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 2351 { 2352 PetscErrorCode ierr; 2353 PetscInt ii; 2354 Mat_SeqSBAIJ *sbaij; 2355 2356 PetscFunctionBegin; 2357 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2358 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2359 2360 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2361 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2362 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2363 ierr = MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 2364 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2365 ierr = PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);CHKERRQ(ierr); 2366 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2367 2368 sbaij->i = i; 2369 sbaij->j = j; 2370 sbaij->a = a; 2371 2372 sbaij->singlemalloc = PETSC_FALSE; 2373 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2374 sbaij->free_a = PETSC_FALSE; 2375 sbaij->free_ij = PETSC_FALSE; 2376 sbaij->free_imax_ilen = PETSC_TRUE; 2377 2378 for (ii=0; ii<m; ii++) { 2379 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2380 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]); 2381 } 2382 if (PetscDefined(USE_DEBUG)) { 2383 for (ii=0; ii<sbaij->i[m]; ii++) { 2384 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2385 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]); 2386 } 2387 } 2388 2389 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2390 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2391 PetscFunctionReturn(0); 2392 } 2393 2394 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2395 { 2396 PetscErrorCode ierr; 2397 PetscMPIInt size; 2398 2399 PetscFunctionBegin; 2400 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 2401 if (size == 1 && scall == MAT_REUSE_MATRIX) { 2402 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2403 } else { 2404 ierr = MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2405 } 2406 PetscFunctionReturn(0); 2407 } 2408