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 591 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 592 { 593 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 594 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 595 PetscInt *ai = a->i,*ailen = a->ilen; 596 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 597 MatScalar *ap,*aa = a->a; 598 599 PetscFunctionBegin; 600 for (k=0; k<m; k++) { /* loop over rows */ 601 row = im[k]; brow = row/bs; 602 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 603 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); 604 rp = aj + ai[brow]; ap = aa + bs2*ai[brow]; 605 nrow = ailen[brow]; 606 for (l=0; l<n; l++) { /* loop over columns */ 607 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 608 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); 609 col = in[l]; 610 bcol = col/bs; 611 cidx = col%bs; 612 ridx = row%bs; 613 high = nrow; 614 low = 0; /* assume unsorted */ 615 while (high-low > 5) { 616 t = (low+high)/2; 617 if (rp[t] > bcol) high = t; 618 else low = t; 619 } 620 for (i=low; i<high; i++) { 621 if (rp[i] > bcol) break; 622 if (rp[i] == bcol) { 623 *v++ = ap[bs2*i+bs*cidx+ridx]; 624 goto finished; 625 } 626 } 627 *v++ = 0.0; 628 finished:; 629 } 630 } 631 PetscFunctionReturn(0); 632 } 633 634 PetscErrorCode MatPermute_SeqSBAIJ(Mat A,IS rowp,IS colp,Mat *B) 635 { 636 Mat C; 637 PetscErrorCode ierr; 638 639 PetscFunctionBegin; 640 ierr = MatConvert(A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 641 ierr = MatPermute(C,rowp,colp,B);CHKERRQ(ierr); 642 ierr = MatDestroy(&C);CHKERRQ(ierr); 643 if (rowp == colp) { 644 ierr = MatConvert(*B,MATSEQSBAIJ,MAT_INPLACE_MATRIX,B);CHKERRQ(ierr); 645 } 646 PetscFunctionReturn(0); 647 } 648 649 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 650 { 651 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 652 PetscErrorCode ierr; 653 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 654 PetscInt *imax =a->imax,*ai=a->i,*ailen=a->ilen; 655 PetscInt *aj =a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 656 PetscBool roworiented=a->roworiented; 657 const PetscScalar *value = v; 658 MatScalar *ap,*aa = a->a,*bap; 659 660 PetscFunctionBegin; 661 if (roworiented) stepval = (n-1)*bs; 662 else stepval = (m-1)*bs; 663 664 for (k=0; k<m; k++) { /* loop over added rows */ 665 row = im[k]; 666 if (row < 0) continue; 667 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); 668 rp = aj + ai[row]; 669 ap = aa + bs2*ai[row]; 670 rmax = imax[row]; 671 nrow = ailen[row]; 672 low = 0; 673 high = nrow; 674 for (l=0; l<n; l++) { /* loop over added columns */ 675 if (in[l] < 0) continue; 676 col = in[l]; 677 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); 678 if (col < row) { 679 if (a->ignore_ltriangular) continue; /* ignore lower triangular block */ 680 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)"); 681 } 682 if (roworiented) value = v + k*(stepval+bs)*bs + l*bs; 683 else value = v + l*(stepval+bs)*bs + k*bs; 684 685 if (col <= lastcol) low = 0; 686 else high = nrow; 687 688 lastcol = col; 689 while (high-low > 7) { 690 t = (low+high)/2; 691 if (rp[t] > col) high = t; 692 else low = t; 693 } 694 for (i=low; i<high; i++) { 695 if (rp[i] > col) break; 696 if (rp[i] == col) { 697 bap = ap + bs2*i; 698 if (roworiented) { 699 if (is == ADD_VALUES) { 700 for (ii=0; ii<bs; ii++,value+=stepval) { 701 for (jj=ii; jj<bs2; jj+=bs) { 702 bap[jj] += *value++; 703 } 704 } 705 } else { 706 for (ii=0; ii<bs; ii++,value+=stepval) { 707 for (jj=ii; jj<bs2; jj+=bs) { 708 bap[jj] = *value++; 709 } 710 } 711 } 712 } else { 713 if (is == ADD_VALUES) { 714 for (ii=0; ii<bs; ii++,value+=stepval) { 715 for (jj=0; jj<bs; jj++) { 716 *bap++ += *value++; 717 } 718 } 719 } else { 720 for (ii=0; ii<bs; ii++,value+=stepval) { 721 for (jj=0; jj<bs; jj++) { 722 *bap++ = *value++; 723 } 724 } 725 } 726 } 727 goto noinsert2; 728 } 729 } 730 if (nonew == 1) goto noinsert2; 731 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); 732 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 733 N = nrow++ - 1; high++; 734 /* shift up all the later entries in this row */ 735 ierr = PetscArraymove(rp+i+1,rp+i,N-i+1);CHKERRQ(ierr); 736 ierr = PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));CHKERRQ(ierr); 737 ierr = PetscArrayzero(ap+bs2*i,bs2);CHKERRQ(ierr); 738 rp[i] = col; 739 bap = ap + bs2*i; 740 if (roworiented) { 741 for (ii=0; ii<bs; ii++,value+=stepval) { 742 for (jj=ii; jj<bs2; jj+=bs) { 743 bap[jj] = *value++; 744 } 745 } 746 } else { 747 for (ii=0; ii<bs; ii++,value+=stepval) { 748 for (jj=0; jj<bs; jj++) { 749 *bap++ = *value++; 750 } 751 } 752 } 753 noinsert2:; 754 low = i; 755 } 756 ailen[row] = nrow; 757 } 758 PetscFunctionReturn(0); 759 } 760 761 /* 762 This is not yet used 763 */ 764 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 765 { 766 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 767 PetscErrorCode ierr; 768 const PetscInt *ai = a->i, *aj = a->j,*cols; 769 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 770 PetscBool flag; 771 772 PetscFunctionBegin; 773 ierr = PetscMalloc1(m,&ns);CHKERRQ(ierr); 774 while (i < m) { 775 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 776 /* Limits the number of elements in a node to 'a->inode.limit' */ 777 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 778 nzy = ai[j+1] - ai[j]; 779 if (nzy != (nzx - j + i)) break; 780 ierr = PetscArraycmp(aj + ai[i] + j - i,aj + ai[j],nzy,&flag);CHKERRQ(ierr); 781 if (!flag) break; 782 } 783 ns[node_count++] = blk_size; 784 785 i = j; 786 } 787 if (!a->inode.size && m && node_count > .9*m) { 788 ierr = PetscFree(ns);CHKERRQ(ierr); 789 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 790 } else { 791 a->inode.node_count = node_count; 792 793 ierr = PetscMalloc1(node_count,&a->inode.size);CHKERRQ(ierr); 794 ierr = PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt));CHKERRQ(ierr); 795 ierr = PetscArraycpy(a->inode.size,ns,node_count);CHKERRQ(ierr); 796 ierr = PetscFree(ns);CHKERRQ(ierr); 797 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 798 799 /* count collections of adjacent columns in each inode */ 800 row = 0; 801 cnt = 0; 802 for (i=0; i<node_count; i++) { 803 cols = aj + ai[row] + a->inode.size[i]; 804 nz = ai[row+1] - ai[row] - a->inode.size[i]; 805 for (j=1; j<nz; j++) { 806 if (cols[j] != cols[j-1]+1) cnt++; 807 } 808 cnt++; 809 row += a->inode.size[i]; 810 } 811 ierr = PetscMalloc1(2*cnt,&counts);CHKERRQ(ierr); 812 cnt = 0; 813 row = 0; 814 for (i=0; i<node_count; i++) { 815 cols = aj + ai[row] + a->inode.size[i]; 816 counts[2*cnt] = cols[0]; 817 nz = ai[row+1] - ai[row] - a->inode.size[i]; 818 cnt2 = 1; 819 for (j=1; j<nz; j++) { 820 if (cols[j] != cols[j-1]+1) { 821 counts[2*(cnt++)+1] = cnt2; 822 counts[2*cnt] = cols[j]; 823 cnt2 = 1; 824 } else cnt2++; 825 } 826 counts[2*(cnt++)+1] = cnt2; 827 row += a->inode.size[i]; 828 } 829 ierr = PetscIntView(2*cnt,counts,NULL);CHKERRQ(ierr); 830 } 831 PetscFunctionReturn(0); 832 } 833 834 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 835 { 836 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 837 PetscErrorCode ierr; 838 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 839 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 840 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 841 MatScalar *aa = a->a,*ap; 842 843 PetscFunctionBegin; 844 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 845 846 if (m) rmax = ailen[0]; 847 for (i=1; i<mbs; i++) { 848 /* move each row back by the amount of empty slots (fshift) before it*/ 849 fshift += imax[i-1] - ailen[i-1]; 850 rmax = PetscMax(rmax,ailen[i]); 851 if (fshift) { 852 ip = aj + ai[i]; 853 ap = aa + bs2*ai[i]; 854 N = ailen[i]; 855 ierr = PetscArraymove(ip-fshift,ip,N);CHKERRQ(ierr); 856 ierr = PetscArraymove(ap-bs2*fshift,ap,bs2*N);CHKERRQ(ierr); 857 } 858 ai[i] = ai[i-1] + ailen[i-1]; 859 } 860 if (mbs) { 861 fshift += imax[mbs-1] - ailen[mbs-1]; 862 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 863 } 864 /* reset ilen and imax for each row */ 865 for (i=0; i<mbs; i++) { 866 ailen[i] = imax[i] = ai[i+1] - ai[i]; 867 } 868 a->nz = ai[mbs]; 869 870 /* diagonals may have moved, reset it */ 871 if (a->diag) { 872 ierr = PetscArraycpy(a->diag,ai,mbs);CHKERRQ(ierr); 873 } 874 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); 875 876 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); 877 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 878 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 879 880 A->info.mallocs += a->reallocs; 881 a->reallocs = 0; 882 A->info.nz_unneeded = (PetscReal)fshift*bs2; 883 a->idiagvalid = PETSC_FALSE; 884 a->rmax = rmax; 885 886 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 887 if (a->jshort && a->free_jshort) { 888 /* when matrix data structure is changed, previous jshort must be replaced */ 889 ierr = PetscFree(a->jshort);CHKERRQ(ierr); 890 } 891 ierr = PetscMalloc1(a->i[A->rmap->n],&a->jshort);CHKERRQ(ierr); 892 ierr = PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short));CHKERRQ(ierr); 893 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 894 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 895 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 896 a->free_jshort = PETSC_TRUE; 897 } 898 PetscFunctionReturn(0); 899 } 900 901 /* 902 This function returns an array of flags which indicate the locations of contiguous 903 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 904 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 905 Assume: sizes should be long enough to hold all the values. 906 */ 907 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 908 { 909 PetscInt i,j,k,row; 910 PetscBool flg; 911 912 PetscFunctionBegin; 913 for (i=0,j=0; i<n; j++) { 914 row = idx[i]; 915 if (row%bs!=0) { /* Not the begining of a block */ 916 sizes[j] = 1; 917 i++; 918 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 919 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 920 i++; 921 } else { /* Begining of the block, so check if the complete block exists */ 922 flg = PETSC_TRUE; 923 for (k=1; k<bs; k++) { 924 if (row+k != idx[i+k]) { /* break in the block */ 925 flg = PETSC_FALSE; 926 break; 927 } 928 } 929 if (flg) { /* No break in the bs */ 930 sizes[j] = bs; 931 i += bs; 932 } else { 933 sizes[j] = 1; 934 i++; 935 } 936 } 937 } 938 *bs_max = j; 939 PetscFunctionReturn(0); 940 } 941 942 943 /* Only add/insert a(i,j) with i<=j (blocks). 944 Any a(i,j) with i>j input by user is ingored. 945 */ 946 947 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 948 { 949 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 950 PetscErrorCode ierr; 951 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 952 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 953 PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 954 PetscInt ridx,cidx,bs2=a->bs2; 955 MatScalar *ap,value,*aa=a->a,*bap; 956 957 PetscFunctionBegin; 958 for (k=0; k<m; k++) { /* loop over added rows */ 959 row = im[k]; /* row number */ 960 brow = row/bs; /* block row number */ 961 if (row < 0) continue; 962 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); 963 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 964 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 965 rmax = imax[brow]; /* maximum space allocated for this row */ 966 nrow = ailen[brow]; /* actual length of this row */ 967 low = 0; 968 high = nrow; 969 for (l=0; l<n; l++) { /* loop over added columns */ 970 if (in[l] < 0) continue; 971 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); 972 col = in[l]; 973 bcol = col/bs; /* block col number */ 974 975 if (brow > bcol) { 976 if (a->ignore_ltriangular) continue; /* ignore lower triangular values */ 977 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)"); 978 } 979 980 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 981 if ((brow==bcol && ridx<=cidx) || (brow<bcol)) { 982 /* element value a(k,l) */ 983 if (roworiented) value = v[l + k*n]; 984 else value = v[k + l*m]; 985 986 /* move pointer bap to a(k,l) quickly and add/insert value */ 987 if (col <= lastcol) low = 0; 988 else high = nrow; 989 990 lastcol = col; 991 while (high-low > 7) { 992 t = (low+high)/2; 993 if (rp[t] > bcol) high = t; 994 else low = t; 995 } 996 for (i=low; i<high; i++) { 997 if (rp[i] > bcol) break; 998 if (rp[i] == bcol) { 999 bap = ap + bs2*i + bs*cidx + ridx; 1000 if (is == ADD_VALUES) *bap += value; 1001 else *bap = value; 1002 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1003 if (brow == bcol && ridx < cidx) { 1004 bap = ap + bs2*i + bs*ridx + cidx; 1005 if (is == ADD_VALUES) *bap += value; 1006 else *bap = value; 1007 } 1008 goto noinsert1; 1009 } 1010 } 1011 1012 if (nonew == 1) goto noinsert1; 1013 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 1014 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 1015 1016 N = nrow++ - 1; high++; 1017 /* shift up all the later entries in this row */ 1018 ierr = PetscArraymove(rp+i+1,rp+i,N-i+1);CHKERRQ(ierr); 1019 ierr = PetscArraymove(ap+bs2*(i+1),ap+bs2*i,bs2*(N-i+1));CHKERRQ(ierr); 1020 ierr = PetscArrayzero(ap+bs2*i,bs2);CHKERRQ(ierr); 1021 rp[i] = bcol; 1022 ap[bs2*i + bs*cidx + ridx] = value; 1023 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1024 if (brow == bcol && ridx < cidx) { 1025 ap[bs2*i + bs*ridx + cidx] = value; 1026 } 1027 A->nonzerostate++; 1028 noinsert1:; 1029 low = i; 1030 } 1031 } /* end of loop over added columns */ 1032 ailen[brow] = nrow; 1033 } /* end of loop over added rows */ 1034 PetscFunctionReturn(0); 1035 } 1036 1037 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1038 { 1039 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1040 Mat outA; 1041 PetscErrorCode ierr; 1042 PetscBool row_identity; 1043 1044 PetscFunctionBegin; 1045 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1046 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1047 if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1048 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()! */ 1049 1050 outA = inA; 1051 inA->factortype = MAT_FACTOR_ICC; 1052 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 1053 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 1054 1055 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1056 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr); 1057 1058 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1059 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1060 a->row = row; 1061 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1062 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1063 a->col = row; 1064 1065 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1066 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1067 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 1068 1069 if (!a->solve_work) { 1070 ierr = PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);CHKERRQ(ierr); 1071 ierr = PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1072 } 1073 1074 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1075 PetscFunctionReturn(0); 1076 } 1077 1078 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1079 { 1080 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data; 1081 PetscInt i,nz,n; 1082 PetscErrorCode ierr; 1083 1084 PetscFunctionBegin; 1085 nz = baij->maxnz; 1086 n = mat->cmap->n; 1087 for (i=0; i<nz; i++) baij->j[i] = indices[i]; 1088 1089 baij->nz = nz; 1090 for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i]; 1091 1092 ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1093 PetscFunctionReturn(0); 1094 } 1095 1096 /*@ 1097 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1098 in the matrix. 1099 1100 Input Parameters: 1101 + mat - the SeqSBAIJ matrix 1102 - indices - the column indices 1103 1104 Level: advanced 1105 1106 Notes: 1107 This can be called if you have precomputed the nonzero structure of the 1108 matrix and want to provide it to the matrix object to improve the performance 1109 of the MatSetValues() operation. 1110 1111 You MUST have set the correct numbers of nonzeros per row in the call to 1112 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1113 1114 MUST be called before any calls to MatSetValues() 1115 1116 .seealso: MatCreateSeqSBAIJ 1117 @*/ 1118 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1119 { 1120 PetscErrorCode ierr; 1121 1122 PetscFunctionBegin; 1123 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1124 PetscValidPointer(indices,2); 1125 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 1126 PetscFunctionReturn(0); 1127 } 1128 1129 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1130 { 1131 PetscErrorCode ierr; 1132 PetscBool isbaij; 1133 1134 PetscFunctionBegin; 1135 ierr = PetscObjectTypeCompareAny((PetscObject)B,&isbaij,MATSEQSBAIJ,MATMPISBAIJ,"");CHKERRQ(ierr); 1136 if (!isbaij) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)B)->type_name); 1137 /* If the two matrices have the same copy implementation and nonzero pattern, use fast copy. */ 1138 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) { 1139 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1140 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1141 1142 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"); 1143 if (a->mbs != b->mbs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of rows in two matrices are different"); 1144 if (a->bs2 != b->bs2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different block size"); 1145 ierr = PetscArraycpy(b->a,a->a,a->bs2*a->i[a->mbs]);CHKERRQ(ierr); 1146 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 1147 } else { 1148 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1149 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1150 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1151 } 1152 PetscFunctionReturn(0); 1153 } 1154 1155 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A) 1156 { 1157 PetscErrorCode ierr; 1158 1159 PetscFunctionBegin; 1160 ierr = MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,NULL);CHKERRQ(ierr); 1161 PetscFunctionReturn(0); 1162 } 1163 1164 static PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1165 { 1166 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1167 1168 PetscFunctionBegin; 1169 *array = a->a; 1170 PetscFunctionReturn(0); 1171 } 1172 1173 static PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1174 { 1175 PetscFunctionBegin; 1176 *array = NULL; 1177 PetscFunctionReturn(0); 1178 } 1179 1180 PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz) 1181 { 1182 PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs; 1183 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data; 1184 Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data; 1185 PetscErrorCode ierr; 1186 1187 PetscFunctionBegin; 1188 /* Set the number of nonzeros in the new matrix */ 1189 ierr = MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 1190 PetscFunctionReturn(0); 1191 } 1192 1193 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1194 { 1195 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data; 1196 PetscErrorCode ierr; 1197 PetscInt bs=Y->rmap->bs,bs2=bs*bs; 1198 PetscBLASInt one = 1; 1199 1200 PetscFunctionBegin; 1201 if (str == SAME_NONZERO_PATTERN) { 1202 PetscScalar alpha = a; 1203 PetscBLASInt bnz; 1204 ierr = PetscBLASIntCast(x->nz*bs2,&bnz);CHKERRQ(ierr); 1205 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 1206 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 1207 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1208 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);CHKERRQ(ierr); 1209 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1210 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);CHKERRQ(ierr); 1211 } else { 1212 Mat B; 1213 PetscInt *nnz; 1214 if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size"); 1215 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1216 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 1217 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 1218 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 1219 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 1220 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 1221 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 1222 ierr = MatSetType(B,((PetscObject)Y)->type_name);CHKERRQ(ierr); 1223 ierr = MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);CHKERRQ(ierr); 1224 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1225 1226 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 1227 1228 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 1229 ierr = PetscFree(nnz);CHKERRQ(ierr); 1230 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1231 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 1232 } 1233 PetscFunctionReturn(0); 1234 } 1235 1236 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1237 { 1238 PetscFunctionBegin; 1239 *flg = PETSC_TRUE; 1240 PetscFunctionReturn(0); 1241 } 1242 1243 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1244 { 1245 PetscFunctionBegin; 1246 *flg = PETSC_TRUE; 1247 PetscFunctionReturn(0); 1248 } 1249 1250 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1251 { 1252 PetscFunctionBegin; 1253 *flg = PETSC_FALSE; 1254 PetscFunctionReturn(0); 1255 } 1256 1257 PetscErrorCode MatConjugate_SeqSBAIJ(Mat A) 1258 { 1259 #if defined(PETSC_USE_COMPLEX) 1260 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1261 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1262 MatScalar *aa = a->a; 1263 1264 PetscFunctionBegin; 1265 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1266 #else 1267 PetscFunctionBegin; 1268 #endif 1269 PetscFunctionReturn(0); 1270 } 1271 1272 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1273 { 1274 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1275 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1276 MatScalar *aa = a->a; 1277 1278 PetscFunctionBegin; 1279 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1280 PetscFunctionReturn(0); 1281 } 1282 1283 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1284 { 1285 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1286 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1287 MatScalar *aa = a->a; 1288 1289 PetscFunctionBegin; 1290 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1291 PetscFunctionReturn(0); 1292 } 1293 1294 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1295 { 1296 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1297 PetscErrorCode ierr; 1298 PetscInt i,j,k,count; 1299 PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col; 1300 PetscScalar zero = 0.0; 1301 MatScalar *aa; 1302 const PetscScalar *xx; 1303 PetscScalar *bb; 1304 PetscBool *zeroed,vecs = PETSC_FALSE; 1305 1306 PetscFunctionBegin; 1307 /* fix right hand side if needed */ 1308 if (x && b) { 1309 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1310 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1311 vecs = PETSC_TRUE; 1312 } 1313 1314 /* zero the columns */ 1315 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1316 for (i=0; i<is_n; i++) { 1317 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]); 1318 zeroed[is_idx[i]] = PETSC_TRUE; 1319 } 1320 if (vecs) { 1321 for (i=0; i<A->rmap->N; i++) { 1322 row = i/bs; 1323 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1324 for (k=0; k<bs; k++) { 1325 col = bs*baij->j[j] + k; 1326 if (col <= i) continue; 1327 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1328 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col]; 1329 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i]; 1330 } 1331 } 1332 } 1333 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]]; 1334 } 1335 1336 for (i=0; i<A->rmap->N; i++) { 1337 if (!zeroed[i]) { 1338 row = i/bs; 1339 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1340 for (k=0; k<bs; k++) { 1341 col = bs*baij->j[j] + k; 1342 if (zeroed[col]) { 1343 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1344 aa[0] = 0.0; 1345 } 1346 } 1347 } 1348 } 1349 } 1350 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1351 if (vecs) { 1352 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1353 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1354 } 1355 1356 /* zero the rows */ 1357 for (i=0; i<is_n; i++) { 1358 row = is_idx[i]; 1359 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1360 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1361 for (k=0; k<count; k++) { 1362 aa[0] = zero; 1363 aa += bs; 1364 } 1365 if (diag != 0.0) { 1366 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1367 } 1368 } 1369 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1370 PetscFunctionReturn(0); 1371 } 1372 1373 PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a) 1374 { 1375 PetscErrorCode ierr; 1376 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data; 1377 1378 PetscFunctionBegin; 1379 if (!Y->preallocated || !aij->nz) { 1380 ierr = MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);CHKERRQ(ierr); 1381 } 1382 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 1383 PetscFunctionReturn(0); 1384 } 1385 1386 /* -------------------------------------------------------------------*/ 1387 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1388 MatGetRow_SeqSBAIJ, 1389 MatRestoreRow_SeqSBAIJ, 1390 MatMult_SeqSBAIJ_N, 1391 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1392 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1393 MatMultAdd_SeqSBAIJ_N, 1394 NULL, 1395 NULL, 1396 NULL, 1397 /* 10*/ NULL, 1398 NULL, 1399 MatCholeskyFactor_SeqSBAIJ, 1400 MatSOR_SeqSBAIJ, 1401 MatTranspose_SeqSBAIJ, 1402 /* 15*/ MatGetInfo_SeqSBAIJ, 1403 MatEqual_SeqSBAIJ, 1404 MatGetDiagonal_SeqSBAIJ, 1405 MatDiagonalScale_SeqSBAIJ, 1406 MatNorm_SeqSBAIJ, 1407 /* 20*/ NULL, 1408 MatAssemblyEnd_SeqSBAIJ, 1409 MatSetOption_SeqSBAIJ, 1410 MatZeroEntries_SeqSBAIJ, 1411 /* 24*/ NULL, 1412 NULL, 1413 NULL, 1414 NULL, 1415 NULL, 1416 /* 29*/ MatSetUp_SeqSBAIJ, 1417 NULL, 1418 NULL, 1419 NULL, 1420 NULL, 1421 /* 34*/ MatDuplicate_SeqSBAIJ, 1422 NULL, 1423 NULL, 1424 NULL, 1425 MatICCFactor_SeqSBAIJ, 1426 /* 39*/ MatAXPY_SeqSBAIJ, 1427 MatCreateSubMatrices_SeqSBAIJ, 1428 MatIncreaseOverlap_SeqSBAIJ, 1429 MatGetValues_SeqSBAIJ, 1430 MatCopy_SeqSBAIJ, 1431 /* 44*/ NULL, 1432 MatScale_SeqSBAIJ, 1433 MatShift_SeqSBAIJ, 1434 NULL, 1435 MatZeroRowsColumns_SeqSBAIJ, 1436 /* 49*/ NULL, 1437 MatGetRowIJ_SeqSBAIJ, 1438 MatRestoreRowIJ_SeqSBAIJ, 1439 NULL, 1440 NULL, 1441 /* 54*/ NULL, 1442 NULL, 1443 NULL, 1444 MatPermute_SeqSBAIJ, 1445 MatSetValuesBlocked_SeqSBAIJ, 1446 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1447 NULL, 1448 NULL, 1449 NULL, 1450 NULL, 1451 /* 64*/ NULL, 1452 NULL, 1453 NULL, 1454 NULL, 1455 NULL, 1456 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1457 NULL, 1458 MatConvert_MPISBAIJ_Basic, 1459 NULL, 1460 NULL, 1461 /* 74*/ NULL, 1462 NULL, 1463 NULL, 1464 NULL, 1465 NULL, 1466 /* 79*/ NULL, 1467 NULL, 1468 NULL, 1469 MatGetInertia_SeqSBAIJ, 1470 MatLoad_SeqSBAIJ, 1471 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1472 MatIsHermitian_SeqSBAIJ, 1473 MatIsStructurallySymmetric_SeqSBAIJ, 1474 NULL, 1475 NULL, 1476 /* 89*/ NULL, 1477 NULL, 1478 NULL, 1479 NULL, 1480 NULL, 1481 /* 94*/ NULL, 1482 NULL, 1483 NULL, 1484 NULL, 1485 NULL, 1486 /* 99*/ NULL, 1487 NULL, 1488 NULL, 1489 MatConjugate_SeqSBAIJ, 1490 NULL, 1491 /*104*/ NULL, 1492 MatRealPart_SeqSBAIJ, 1493 MatImaginaryPart_SeqSBAIJ, 1494 MatGetRowUpperTriangular_SeqSBAIJ, 1495 MatRestoreRowUpperTriangular_SeqSBAIJ, 1496 /*109*/ NULL, 1497 NULL, 1498 NULL, 1499 NULL, 1500 MatMissingDiagonal_SeqSBAIJ, 1501 /*114*/ NULL, 1502 NULL, 1503 NULL, 1504 NULL, 1505 NULL, 1506 /*119*/ NULL, 1507 NULL, 1508 NULL, 1509 NULL, 1510 NULL, 1511 /*124*/ NULL, 1512 NULL, 1513 NULL, 1514 NULL, 1515 NULL, 1516 /*129*/ NULL, 1517 NULL, 1518 NULL, 1519 NULL, 1520 NULL, 1521 /*134*/ NULL, 1522 NULL, 1523 NULL, 1524 NULL, 1525 NULL, 1526 /*139*/ MatSetBlockSizes_Default, 1527 NULL, 1528 NULL, 1529 NULL, 1530 NULL, 1531 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ 1532 }; 1533 1534 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1535 { 1536 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1537 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1538 PetscErrorCode ierr; 1539 1540 PetscFunctionBegin; 1541 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1542 1543 /* allocate space for values if not already there */ 1544 if (!aij->saved_values) { 1545 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 1546 } 1547 1548 /* copy values over */ 1549 ierr = PetscArraycpy(aij->saved_values,aij->a,nz);CHKERRQ(ierr); 1550 PetscFunctionReturn(0); 1551 } 1552 1553 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1554 { 1555 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1556 PetscErrorCode ierr; 1557 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1558 1559 PetscFunctionBegin; 1560 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1561 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1562 1563 /* copy values over */ 1564 ierr = PetscArraycpy(aij->a,aij->saved_values,nz);CHKERRQ(ierr); 1565 PetscFunctionReturn(0); 1566 } 1567 1568 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1569 { 1570 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1571 PetscErrorCode ierr; 1572 PetscInt i,mbs,nbs,bs2; 1573 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1574 1575 PetscFunctionBegin; 1576 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1577 1578 ierr = MatSetBlockSize(B,PetscAbs(bs));CHKERRQ(ierr); 1579 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1580 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1581 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); 1582 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1583 1584 B->preallocated = PETSC_TRUE; 1585 1586 mbs = B->rmap->N/bs; 1587 nbs = B->cmap->n/bs; 1588 bs2 = bs*bs; 1589 1590 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"); 1591 1592 if (nz == MAT_SKIP_ALLOCATION) { 1593 skipallocation = PETSC_TRUE; 1594 nz = 0; 1595 } 1596 1597 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1598 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1599 if (nnz) { 1600 for (i=0; i<mbs; i++) { 1601 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]); 1602 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); 1603 } 1604 } 1605 1606 B->ops->mult = MatMult_SeqSBAIJ_N; 1607 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1608 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1609 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1610 1611 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1612 if (!flg) { 1613 switch (bs) { 1614 case 1: 1615 B->ops->mult = MatMult_SeqSBAIJ_1; 1616 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1617 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1618 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1619 break; 1620 case 2: 1621 B->ops->mult = MatMult_SeqSBAIJ_2; 1622 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1623 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1624 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1625 break; 1626 case 3: 1627 B->ops->mult = MatMult_SeqSBAIJ_3; 1628 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1629 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1630 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1631 break; 1632 case 4: 1633 B->ops->mult = MatMult_SeqSBAIJ_4; 1634 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1635 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1636 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1637 break; 1638 case 5: 1639 B->ops->mult = MatMult_SeqSBAIJ_5; 1640 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1641 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1642 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1643 break; 1644 case 6: 1645 B->ops->mult = MatMult_SeqSBAIJ_6; 1646 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1647 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1648 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1649 break; 1650 case 7: 1651 B->ops->mult = MatMult_SeqSBAIJ_7; 1652 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1653 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1654 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1655 break; 1656 } 1657 } 1658 1659 b->mbs = mbs; 1660 b->nbs = nbs; 1661 if (!skipallocation) { 1662 if (!b->imax) { 1663 ierr = PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);CHKERRQ(ierr); 1664 1665 b->free_imax_ilen = PETSC_TRUE; 1666 1667 ierr = PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1668 } 1669 if (!nnz) { 1670 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1671 else if (nz <= 0) nz = 1; 1672 nz = PetscMin(nbs,nz); 1673 for (i=0; i<mbs; i++) b->imax[i] = nz; 1674 ierr = PetscIntMultError(nz,mbs,&nz);CHKERRQ(ierr); 1675 } else { 1676 PetscInt64 nz64 = 0; 1677 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 1678 ierr = PetscIntCast(nz64,&nz);CHKERRQ(ierr); 1679 } 1680 /* b->ilen will count nonzeros in each block row so far. */ 1681 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1682 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1683 1684 /* allocate the matrix space */ 1685 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1686 ierr = PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);CHKERRQ(ierr); 1687 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1688 ierr = PetscArrayzero(b->a,nz*bs2);CHKERRQ(ierr); 1689 ierr = PetscArrayzero(b->j,nz);CHKERRQ(ierr); 1690 1691 b->singlemalloc = PETSC_TRUE; 1692 1693 /* pointer to beginning of each row */ 1694 b->i[0] = 0; 1695 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1696 1697 b->free_a = PETSC_TRUE; 1698 b->free_ij = PETSC_TRUE; 1699 } else { 1700 b->free_a = PETSC_FALSE; 1701 b->free_ij = PETSC_FALSE; 1702 } 1703 1704 b->bs2 = bs2; 1705 b->nz = 0; 1706 b->maxnz = nz; 1707 b->inew = NULL; 1708 b->jnew = NULL; 1709 b->anew = NULL; 1710 b->a2anew = NULL; 1711 b->permute = PETSC_FALSE; 1712 1713 B->was_assembled = PETSC_FALSE; 1714 B->assembled = PETSC_FALSE; 1715 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1716 PetscFunctionReturn(0); 1717 } 1718 1719 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1720 { 1721 PetscInt i,j,m,nz,anz, nz_max=0,*nnz; 1722 PetscScalar *values=NULL; 1723 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1724 PetscErrorCode ierr; 1725 1726 PetscFunctionBegin; 1727 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs); 1728 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1729 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1730 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1731 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1732 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1733 m = B->rmap->n/bs; 1734 1735 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]); 1736 ierr = PetscMalloc1(m+1,&nnz);CHKERRQ(ierr); 1737 for (i=0; i<m; i++) { 1738 nz = ii[i+1] - ii[i]; 1739 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz); 1740 anz = 0; 1741 for (j=0; j<nz; j++) { 1742 /* count only values on the diagonal or above */ 1743 if (jj[ii[i] + j] >= i) { 1744 anz = nz - j; 1745 break; 1746 } 1747 } 1748 nz_max = PetscMax(nz_max,anz); 1749 nnz[i] = anz; 1750 } 1751 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1752 ierr = PetscFree(nnz);CHKERRQ(ierr); 1753 1754 values = (PetscScalar*)V; 1755 if (!values) { 1756 ierr = PetscCalloc1(bs*bs*nz_max,&values);CHKERRQ(ierr); 1757 } 1758 for (i=0; i<m; i++) { 1759 PetscInt ncols = ii[i+1] - ii[i]; 1760 const PetscInt *icols = jj + ii[i]; 1761 if (!roworiented || bs == 1) { 1762 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1763 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr); 1764 } else { 1765 for (j=0; j<ncols; j++) { 1766 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1767 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr); 1768 } 1769 } 1770 } 1771 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); } 1772 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1773 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1774 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1775 PetscFunctionReturn(0); 1776 } 1777 1778 /* 1779 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1780 */ 1781 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1782 { 1783 PetscErrorCode ierr; 1784 PetscBool flg = PETSC_FALSE; 1785 PetscInt bs = B->rmap->bs; 1786 1787 PetscFunctionBegin; 1788 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1789 if (flg) bs = 8; 1790 1791 if (!natural) { 1792 switch (bs) { 1793 case 1: 1794 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1795 break; 1796 case 2: 1797 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1798 break; 1799 case 3: 1800 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1801 break; 1802 case 4: 1803 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1804 break; 1805 case 5: 1806 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1807 break; 1808 case 6: 1809 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1810 break; 1811 case 7: 1812 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1813 break; 1814 default: 1815 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1816 break; 1817 } 1818 } else { 1819 switch (bs) { 1820 case 1: 1821 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1822 break; 1823 case 2: 1824 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1825 break; 1826 case 3: 1827 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1828 break; 1829 case 4: 1830 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1831 break; 1832 case 5: 1833 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1834 break; 1835 case 6: 1836 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1837 break; 1838 case 7: 1839 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1840 break; 1841 default: 1842 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1843 break; 1844 } 1845 } 1846 PetscFunctionReturn(0); 1847 } 1848 1849 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat,MatType,MatReuse,Mat*); 1850 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat,MatType,MatReuse,Mat*); 1851 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type) 1852 { 1853 PetscFunctionBegin; 1854 *type = MATSOLVERPETSC; 1855 PetscFunctionReturn(0); 1856 } 1857 1858 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1859 { 1860 PetscInt n = A->rmap->n; 1861 PetscErrorCode ierr; 1862 1863 PetscFunctionBegin; 1864 #if defined(PETSC_USE_COMPLEX) 1865 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"); 1866 #endif 1867 1868 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 1869 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1870 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1871 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1872 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1873 1874 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1875 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1876 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);CHKERRQ(ierr); 1877 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);CHKERRQ(ierr); 1878 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1879 1880 (*B)->factortype = ftype; 1881 (*B)->canuseordering = PETSC_TRUE; 1882 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 1883 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 1884 ierr = PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc);CHKERRQ(ierr); 1885 PetscFunctionReturn(0); 1886 } 1887 1888 /*@C 1889 MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored 1890 1891 Not Collective 1892 1893 Input Parameter: 1894 . mat - a MATSEQSBAIJ matrix 1895 1896 Output Parameter: 1897 . array - pointer to the data 1898 1899 Level: intermediate 1900 1901 .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray() 1902 @*/ 1903 PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array) 1904 { 1905 PetscErrorCode ierr; 1906 1907 PetscFunctionBegin; 1908 ierr = PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1909 PetscFunctionReturn(0); 1910 } 1911 1912 /*@C 1913 MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray() 1914 1915 Not Collective 1916 1917 Input Parameters: 1918 + mat - a MATSEQSBAIJ matrix 1919 - array - pointer to the data 1920 1921 Level: intermediate 1922 1923 .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray() 1924 @*/ 1925 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array) 1926 { 1927 PetscErrorCode ierr; 1928 1929 PetscFunctionBegin; 1930 ierr = PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1931 PetscFunctionReturn(0); 1932 } 1933 1934 /*MC 1935 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1936 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1937 1938 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1939 can call MatSetOption(Mat, MAT_HERMITIAN). 1940 1941 Options Database Keys: 1942 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1943 1944 Notes: 1945 By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1946 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1947 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. 1948 1949 The number of rows in the matrix must be less than or equal to the number of columns 1950 1951 Level: beginner 1952 1953 .seealso: MatCreateSeqSBAIJ(), MatType, MATMPISBAIJ 1954 M*/ 1955 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1956 { 1957 Mat_SeqSBAIJ *b; 1958 PetscErrorCode ierr; 1959 PetscMPIInt size; 1960 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1961 1962 PetscFunctionBegin; 1963 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRMPI(ierr); 1964 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1965 1966 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1967 B->data = (void*)b; 1968 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1969 1970 B->ops->destroy = MatDestroy_SeqSBAIJ; 1971 B->ops->view = MatView_SeqSBAIJ; 1972 b->row = NULL; 1973 b->icol = NULL; 1974 b->reallocs = 0; 1975 b->saved_values = NULL; 1976 b->inode.limit = 5; 1977 b->inode.max_limit = 5; 1978 1979 b->roworiented = PETSC_TRUE; 1980 b->nonew = 0; 1981 b->diag = NULL; 1982 b->solve_work = NULL; 1983 b->mult_work = NULL; 1984 B->spptr = NULL; 1985 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1986 b->keepnonzeropattern = PETSC_FALSE; 1987 1988 b->inew = NULL; 1989 b->jnew = NULL; 1990 b->anew = NULL; 1991 b->a2anew = NULL; 1992 b->permute = PETSC_FALSE; 1993 1994 b->ignore_ltriangular = PETSC_TRUE; 1995 1996 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1997 1998 b->getrow_utriangular = PETSC_FALSE; 1999 2000 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 2001 2002 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);CHKERRQ(ierr); 2003 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);CHKERRQ(ierr); 2004 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 2005 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 2006 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 2007 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 2008 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 2009 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 2010 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr); 2011 #if defined(PETSC_HAVE_ELEMENTAL) 2012 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);CHKERRQ(ierr); 2013 #endif 2014 #if defined(PETSC_HAVE_SCALAPACK) 2015 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_scalapack_C",MatConvert_SBAIJ_ScaLAPACK);CHKERRQ(ierr); 2016 #endif 2017 2018 B->symmetric = PETSC_TRUE; 2019 B->structurally_symmetric = PETSC_TRUE; 2020 B->symmetric_set = PETSC_TRUE; 2021 B->structurally_symmetric_set = PETSC_TRUE; 2022 B->symmetric_eternal = PETSC_TRUE; 2023 #if defined(PETSC_USE_COMPLEX) 2024 B->hermitian = PETSC_FALSE; 2025 B->hermitian_set = PETSC_FALSE; 2026 #else 2027 B->hermitian = PETSC_TRUE; 2028 B->hermitian_set = PETSC_TRUE; 2029 #endif 2030 2031 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 2032 2033 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 2034 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 2035 if (no_unroll) { 2036 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 2037 } 2038 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 2039 if (no_inode) { 2040 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 2041 } 2042 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 2043 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2044 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 2045 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 2046 PetscFunctionReturn(0); 2047 } 2048 2049 /*@C 2050 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2051 compressed row) format. For good matrix assembly performance the 2052 user should preallocate the matrix storage by setting the parameter nz 2053 (or the array nnz). By setting these parameters accurately, performance 2054 during matrix assembly can be increased by more than a factor of 50. 2055 2056 Collective on Mat 2057 2058 Input Parameters: 2059 + B - the symmetric matrix 2060 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2061 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2062 . nz - number of block nonzeros per block row (same for all rows) 2063 - nnz - array containing the number of block nonzeros in the upper triangular plus 2064 diagonal portion of each block (possibly different for each block row) or NULL 2065 2066 Options Database Keys: 2067 + -mat_no_unroll - uses code that does not unroll the loops in the 2068 block calculations (much slower) 2069 - -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2070 2071 Level: intermediate 2072 2073 Notes: 2074 Specify the preallocated storage with either nz or nnz (not both). 2075 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2076 allocation. See Users-Manual: ch_mat for details. 2077 2078 You can call MatGetInfo() to get information on how effective the preallocation was; 2079 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2080 You can also run with the option -info and look for messages with the string 2081 malloc in them to see if additional memory allocation was needed. 2082 2083 If the nnz parameter is given then the nz parameter is ignored 2084 2085 2086 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2087 @*/ 2088 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2089 { 2090 PetscErrorCode ierr; 2091 2092 PetscFunctionBegin; 2093 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2094 PetscValidType(B,1); 2095 PetscValidLogicalCollectiveInt(B,bs,2); 2096 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2097 PetscFunctionReturn(0); 2098 } 2099 2100 /*@C 2101 MatSeqSBAIJSetPreallocationCSR - Creates a sparse parallel matrix in SBAIJ format using the given nonzero structure and (optional) numerical values 2102 2103 Input Parameters: 2104 + B - the matrix 2105 . bs - size of block, the blocks are ALWAYS square. 2106 . i - the indices into j for the start of each local row (starts with zero) 2107 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2108 - v - optional values in the matrix 2109 2110 Level: advanced 2111 2112 Notes: 2113 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2114 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2115 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2116 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2117 block column and the second index is over columns within a block. 2118 2119 Any entries below the diagonal are ignored 2120 2121 Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries 2122 and usually the numerical values as well 2123 2124 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ 2125 @*/ 2126 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2127 { 2128 PetscErrorCode ierr; 2129 2130 PetscFunctionBegin; 2131 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2132 PetscValidType(B,1); 2133 PetscValidLogicalCollectiveInt(B,bs,2); 2134 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2135 PetscFunctionReturn(0); 2136 } 2137 2138 /*@C 2139 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2140 compressed row) format. For good matrix assembly performance the 2141 user should preallocate the matrix storage by setting the parameter nz 2142 (or the array nnz). By setting these parameters accurately, performance 2143 during matrix assembly can be increased by more than a factor of 50. 2144 2145 Collective 2146 2147 Input Parameters: 2148 + comm - MPI communicator, set to PETSC_COMM_SELF 2149 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2150 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2151 . m - number of rows, or number of columns 2152 . nz - number of block nonzeros per block row (same for all rows) 2153 - nnz - array containing the number of block nonzeros in the upper triangular plus 2154 diagonal portion of each block (possibly different for each block row) or NULL 2155 2156 Output Parameter: 2157 . A - the symmetric matrix 2158 2159 Options Database Keys: 2160 + -mat_no_unroll - uses code that does not unroll the loops in the 2161 block calculations (much slower) 2162 - -mat_block_size - size of the blocks to use 2163 2164 Level: intermediate 2165 2166 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2167 MatXXXXSetPreallocation() paradigm instead of this routine directly. 2168 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2169 2170 Notes: 2171 The number of rows and columns must be divisible by blocksize. 2172 This matrix type does not support complex Hermitian operation. 2173 2174 Specify the preallocated storage with either nz or nnz (not both). 2175 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2176 allocation. See Users-Manual: ch_mat for details. 2177 2178 If the nnz parameter is given then the nz parameter is ignored 2179 2180 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2181 @*/ 2182 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2183 { 2184 PetscErrorCode ierr; 2185 2186 PetscFunctionBegin; 2187 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2188 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2189 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2190 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2191 PetscFunctionReturn(0); 2192 } 2193 2194 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2195 { 2196 Mat C; 2197 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2198 PetscErrorCode ierr; 2199 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2200 2201 PetscFunctionBegin; 2202 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2203 2204 *B = NULL; 2205 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2206 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2207 ierr = MatSetBlockSizesFromMats(C,A,A);CHKERRQ(ierr); 2208 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2209 c = (Mat_SeqSBAIJ*)C->data; 2210 2211 C->preallocated = PETSC_TRUE; 2212 C->factortype = A->factortype; 2213 c->row = NULL; 2214 c->icol = NULL; 2215 c->saved_values = NULL; 2216 c->keepnonzeropattern = a->keepnonzeropattern; 2217 C->assembled = PETSC_TRUE; 2218 2219 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2220 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2221 c->bs2 = a->bs2; 2222 c->mbs = a->mbs; 2223 c->nbs = a->nbs; 2224 2225 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2226 c->imax = a->imax; 2227 c->ilen = a->ilen; 2228 c->free_imax_ilen = PETSC_FALSE; 2229 } else { 2230 ierr = PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);CHKERRQ(ierr); 2231 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2232 for (i=0; i<mbs; i++) { 2233 c->imax[i] = a->imax[i]; 2234 c->ilen[i] = a->ilen[i]; 2235 } 2236 c->free_imax_ilen = PETSC_TRUE; 2237 } 2238 2239 /* allocate the matrix space */ 2240 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2241 ierr = PetscMalloc1(bs2*nz,&c->a);CHKERRQ(ierr); 2242 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2243 c->i = a->i; 2244 c->j = a->j; 2245 c->singlemalloc = PETSC_FALSE; 2246 c->free_a = PETSC_TRUE; 2247 c->free_ij = PETSC_FALSE; 2248 c->parent = A; 2249 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2250 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2251 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2252 } else { 2253 ierr = PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);CHKERRQ(ierr); 2254 ierr = PetscArraycpy(c->i,a->i,mbs+1);CHKERRQ(ierr); 2255 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2256 c->singlemalloc = PETSC_TRUE; 2257 c->free_a = PETSC_TRUE; 2258 c->free_ij = PETSC_TRUE; 2259 } 2260 if (mbs > 0) { 2261 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2262 ierr = PetscArraycpy(c->j,a->j,nz);CHKERRQ(ierr); 2263 } 2264 if (cpvalues == MAT_COPY_VALUES) { 2265 ierr = PetscArraycpy(c->a,a->a,bs2*nz);CHKERRQ(ierr); 2266 } else { 2267 ierr = PetscArrayzero(c->a,bs2*nz);CHKERRQ(ierr); 2268 } 2269 if (a->jshort) { 2270 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2271 /* if the parent matrix is reassembled, this child matrix will never notice */ 2272 ierr = PetscMalloc1(nz,&c->jshort);CHKERRQ(ierr); 2273 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2274 ierr = PetscArraycpy(c->jshort,a->jshort,nz);CHKERRQ(ierr); 2275 2276 c->free_jshort = PETSC_TRUE; 2277 } 2278 } 2279 2280 c->roworiented = a->roworiented; 2281 c->nonew = a->nonew; 2282 2283 if (a->diag) { 2284 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2285 c->diag = a->diag; 2286 c->free_diag = PETSC_FALSE; 2287 } else { 2288 ierr = PetscMalloc1(mbs,&c->diag);CHKERRQ(ierr); 2289 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2290 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2291 c->free_diag = PETSC_TRUE; 2292 } 2293 } 2294 c->nz = a->nz; 2295 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2296 c->solve_work = NULL; 2297 c->mult_work = NULL; 2298 2299 *B = C; 2300 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2301 PetscFunctionReturn(0); 2302 } 2303 2304 /* Used for both SeqBAIJ and SeqSBAIJ matrices */ 2305 #define MatLoad_SeqSBAIJ_Binary MatLoad_SeqBAIJ_Binary 2306 2307 PetscErrorCode MatLoad_SeqSBAIJ(Mat mat,PetscViewer viewer) 2308 { 2309 PetscErrorCode ierr; 2310 PetscBool isbinary; 2311 2312 PetscFunctionBegin; 2313 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 2314 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); 2315 ierr = MatLoad_SeqSBAIJ_Binary(mat,viewer);CHKERRQ(ierr); 2316 PetscFunctionReturn(0); 2317 } 2318 2319 /*@ 2320 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2321 (upper triangular entries in CSR format) provided by the user. 2322 2323 Collective 2324 2325 Input Parameters: 2326 + comm - must be an MPI communicator of size 1 2327 . bs - size of block 2328 . m - number of rows 2329 . n - number of columns 2330 . 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 2331 . j - column indices 2332 - a - matrix values 2333 2334 Output Parameter: 2335 . mat - the matrix 2336 2337 Level: advanced 2338 2339 Notes: 2340 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2341 once the matrix is destroyed 2342 2343 You cannot set new nonzero locations into this matrix, that will generate an error. 2344 2345 The i and j indices are 0 based 2346 2347 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 2348 it is the regular CSR format excluding the lower triangular elements. 2349 2350 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2351 2352 @*/ 2353 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 2354 { 2355 PetscErrorCode ierr; 2356 PetscInt ii; 2357 Mat_SeqSBAIJ *sbaij; 2358 2359 PetscFunctionBegin; 2360 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2361 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2362 2363 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2364 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2365 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2366 ierr = MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 2367 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2368 ierr = PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);CHKERRQ(ierr); 2369 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2370 2371 sbaij->i = i; 2372 sbaij->j = j; 2373 sbaij->a = a; 2374 2375 sbaij->singlemalloc = PETSC_FALSE; 2376 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2377 sbaij->free_a = PETSC_FALSE; 2378 sbaij->free_ij = PETSC_FALSE; 2379 sbaij->free_imax_ilen = PETSC_TRUE; 2380 2381 for (ii=0; ii<m; ii++) { 2382 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2383 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]); 2384 } 2385 if (PetscDefined(USE_DEBUG)) { 2386 for (ii=0; ii<sbaij->i[m]; ii++) { 2387 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2388 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]); 2389 } 2390 } 2391 2392 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2393 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2394 PetscFunctionReturn(0); 2395 } 2396 2397 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2398 { 2399 PetscErrorCode ierr; 2400 PetscMPIInt size; 2401 2402 PetscFunctionBegin; 2403 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 2404 if (size == 1 && scall == MAT_REUSE_MATRIX) { 2405 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2406 } else { 2407 ierr = MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2408 } 2409 PetscFunctionReturn(0); 2410 } 2411