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