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