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