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 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqsbstrm_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 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 238 break; 239 case MAT_HERMITIAN: 240 #if defined(PETSC_USE_COMPLEX) /* MAT_HERMITIAN is a synonym for MAT_SYMMETRIC with reals */ 241 if (!A->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call MatAssemblyEnd() first"); 242 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 243 A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian_ushort; 244 } else if (A->cmap->bs == 1) { 245 A->ops->mult = MatMult_SeqSBAIJ_1_Hermitian; 246 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for Hermitian with block size greater than 1"); 247 #endif 248 break; 249 case MAT_SPD: 250 /* These options are handled directly by MatSetOption() */ 251 break; 252 case MAT_SYMMETRIC: 253 case MAT_STRUCTURALLY_SYMMETRIC: 254 case MAT_SYMMETRY_ETERNAL: 255 case MAT_STRUCTURE_ONLY: 256 /* These options are handled directly by MatSetOption() */ 257 break; 258 case MAT_IGNORE_LOWER_TRIANGULAR: 259 a->ignore_ltriangular = flg; 260 break; 261 case MAT_ERROR_LOWER_TRIANGULAR: 262 a->ignore_ltriangular = flg; 263 break; 264 case MAT_GETROW_UPPERTRIANGULAR: 265 a->getrow_utriangular = flg; 266 break; 267 case MAT_SUBMAT_SINGLEIS: 268 break; 269 default: 270 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 271 } 272 PetscFunctionReturn(0); 273 } 274 275 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 276 { 277 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 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()"); 282 283 /* Get the upper triangular part of the row */ 284 ierr = MatGetRow_SeqBAIJ_private(A,row,nz,idx,v,a->i,a->j,a->a);CHKERRQ(ierr); 285 PetscFunctionReturn(0); 286 } 287 288 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 289 { 290 PetscErrorCode ierr; 291 292 PetscFunctionBegin; 293 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);} 294 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);} 295 PetscFunctionReturn(0); 296 } 297 298 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 299 { 300 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 301 302 PetscFunctionBegin; 303 a->getrow_utriangular = PETSC_TRUE; 304 PetscFunctionReturn(0); 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 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->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->rmap->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 high = nrow; 966 lastcol = col; 967 while (high-low > 7) { 968 t = (low+high)/2; 969 if (rp[t] > bcol) high = t; 970 else low = t; 971 } 972 for (i=low; i<high; i++) { 973 if (rp[i] > bcol) break; 974 if (rp[i] == bcol) { 975 bap = ap + bs2*i + bs*cidx + ridx; 976 if (is == ADD_VALUES) *bap += value; 977 else *bap = value; 978 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 979 if (brow == bcol && ridx < cidx) { 980 bap = ap + bs2*i + bs*ridx + cidx; 981 if (is == ADD_VALUES) *bap += value; 982 else *bap = value; 983 } 984 goto noinsert1; 985 } 986 } 987 988 if (nonew == 1) goto noinsert1; 989 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 990 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 991 992 N = nrow++ - 1; high++; 993 /* shift up all the later entries in this row */ 994 for (ii=N; ii>=i; ii--) { 995 rp[ii+1] = rp[ii]; 996 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 997 } 998 if (N>=i) { 999 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1000 } 1001 rp[i] = bcol; 1002 ap[bs2*i + bs*cidx + ridx] = value; 1003 A->nonzerostate++; 1004 noinsert1:; 1005 low = i; 1006 } 1007 } /* end of loop over added columns */ 1008 ailen[brow] = nrow; 1009 } /* end of loop over added rows */ 1010 PetscFunctionReturn(0); 1011 } 1012 1013 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 1014 { 1015 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 1016 Mat outA; 1017 PetscErrorCode ierr; 1018 PetscBool row_identity; 1019 1020 PetscFunctionBegin; 1021 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1022 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1023 if (!row_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1024 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()! */ 1025 1026 outA = inA; 1027 inA->factortype = MAT_FACTOR_ICC; 1028 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 1029 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 1030 1031 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1032 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr); 1033 1034 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1035 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1036 a->row = row; 1037 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1038 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1039 a->col = row; 1040 1041 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1042 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1043 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 1044 1045 if (!a->solve_work) { 1046 ierr = PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);CHKERRQ(ierr); 1047 ierr = PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1048 } 1049 1050 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1051 PetscFunctionReturn(0); 1052 } 1053 1054 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1055 { 1056 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data; 1057 PetscInt i,nz,n; 1058 PetscErrorCode ierr; 1059 1060 PetscFunctionBegin; 1061 nz = baij->maxnz; 1062 n = mat->cmap->n; 1063 for (i=0; i<nz; i++) baij->j[i] = indices[i]; 1064 1065 baij->nz = nz; 1066 for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i]; 1067 1068 ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1069 PetscFunctionReturn(0); 1070 } 1071 1072 /*@ 1073 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1074 in the matrix. 1075 1076 Input Parameters: 1077 + mat - the SeqSBAIJ matrix 1078 - indices - the column indices 1079 1080 Level: advanced 1081 1082 Notes: 1083 This can be called if you have precomputed the nonzero structure of the 1084 matrix and want to provide it to the matrix object to improve the performance 1085 of the MatSetValues() operation. 1086 1087 You MUST have set the correct numbers of nonzeros per row in the call to 1088 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1089 1090 MUST be called before any calls to MatSetValues() 1091 1092 .seealso: MatCreateSeqSBAIJ 1093 @*/ 1094 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1095 { 1096 PetscErrorCode ierr; 1097 1098 PetscFunctionBegin; 1099 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1100 PetscValidPointer(indices,2); 1101 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 1102 PetscFunctionReturn(0); 1103 } 1104 1105 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1106 { 1107 PetscErrorCode ierr; 1108 1109 PetscFunctionBegin; 1110 /* If the two matrices have the same copy implementation, use fast copy. */ 1111 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1112 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1113 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1114 1115 if (a->i[A->rmap->N] != b->i[B->rmap->N]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1116 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1117 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 1118 } else { 1119 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1120 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1121 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1122 } 1123 PetscFunctionReturn(0); 1124 } 1125 1126 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A) 1127 { 1128 PetscErrorCode ierr; 1129 1130 PetscFunctionBegin; 1131 ierr = MatSeqSBAIJSetPreallocation(A,A->rmap->bs,PETSC_DEFAULT,0);CHKERRQ(ierr); 1132 PetscFunctionReturn(0); 1133 } 1134 1135 PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1136 { 1137 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1138 1139 PetscFunctionBegin; 1140 *array = a->a; 1141 PetscFunctionReturn(0); 1142 } 1143 1144 PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1145 { 1146 PetscFunctionBegin; 1147 PetscFunctionReturn(0); 1148 } 1149 1150 PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz) 1151 { 1152 PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs; 1153 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data; 1154 Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data; 1155 PetscErrorCode ierr; 1156 1157 PetscFunctionBegin; 1158 /* Set the number of nonzeros in the new matrix */ 1159 ierr = MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 1160 PetscFunctionReturn(0); 1161 } 1162 1163 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1164 { 1165 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data; 1166 PetscErrorCode ierr; 1167 PetscInt bs=Y->rmap->bs,bs2=bs*bs; 1168 PetscBLASInt one = 1; 1169 1170 PetscFunctionBegin; 1171 if (str == SAME_NONZERO_PATTERN) { 1172 PetscScalar alpha = a; 1173 PetscBLASInt bnz; 1174 ierr = PetscBLASIntCast(x->nz*bs2,&bnz);CHKERRQ(ierr); 1175 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 1176 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 1177 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1178 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);CHKERRQ(ierr); 1179 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1180 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);CHKERRQ(ierr); 1181 } else { 1182 Mat B; 1183 PetscInt *nnz; 1184 if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size"); 1185 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1186 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 1187 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 1188 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 1189 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 1190 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 1191 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 1192 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 1193 ierr = MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);CHKERRQ(ierr); 1194 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1195 1196 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 1197 1198 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 1199 ierr = PetscFree(nnz);CHKERRQ(ierr); 1200 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1201 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 1202 } 1203 PetscFunctionReturn(0); 1204 } 1205 1206 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1207 { 1208 PetscFunctionBegin; 1209 *flg = PETSC_TRUE; 1210 PetscFunctionReturn(0); 1211 } 1212 1213 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1214 { 1215 PetscFunctionBegin; 1216 *flg = PETSC_TRUE; 1217 PetscFunctionReturn(0); 1218 } 1219 1220 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1221 { 1222 PetscFunctionBegin; 1223 *flg = PETSC_FALSE; 1224 PetscFunctionReturn(0); 1225 } 1226 1227 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1228 { 1229 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1230 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1231 MatScalar *aa = a->a; 1232 1233 PetscFunctionBegin; 1234 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1235 PetscFunctionReturn(0); 1236 } 1237 1238 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1239 { 1240 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1241 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1242 MatScalar *aa = a->a; 1243 1244 PetscFunctionBegin; 1245 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1246 PetscFunctionReturn(0); 1247 } 1248 1249 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1250 { 1251 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1252 PetscErrorCode ierr; 1253 PetscInt i,j,k,count; 1254 PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col; 1255 PetscScalar zero = 0.0; 1256 MatScalar *aa; 1257 const PetscScalar *xx; 1258 PetscScalar *bb; 1259 PetscBool *zeroed,vecs = PETSC_FALSE; 1260 1261 PetscFunctionBegin; 1262 /* fix right hand side if needed */ 1263 if (x && b) { 1264 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1265 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1266 vecs = PETSC_TRUE; 1267 } 1268 1269 /* zero the columns */ 1270 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1271 for (i=0; i<is_n; i++) { 1272 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]); 1273 zeroed[is_idx[i]] = PETSC_TRUE; 1274 } 1275 if (vecs) { 1276 for (i=0; i<A->rmap->N; i++) { 1277 row = i/bs; 1278 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1279 for (k=0; k<bs; k++) { 1280 col = bs*baij->j[j] + k; 1281 if (col <= i) continue; 1282 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1283 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col]; 1284 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i]; 1285 } 1286 } 1287 } 1288 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]]; 1289 } 1290 1291 for (i=0; i<A->rmap->N; i++) { 1292 if (!zeroed[i]) { 1293 row = i/bs; 1294 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1295 for (k=0; k<bs; k++) { 1296 col = bs*baij->j[j] + k; 1297 if (zeroed[col]) { 1298 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1299 aa[0] = 0.0; 1300 } 1301 } 1302 } 1303 } 1304 } 1305 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1306 if (vecs) { 1307 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1308 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1309 } 1310 1311 /* zero the rows */ 1312 for (i=0; i<is_n; i++) { 1313 row = is_idx[i]; 1314 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1315 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1316 for (k=0; k<count; k++) { 1317 aa[0] = zero; 1318 aa += bs; 1319 } 1320 if (diag != 0.0) { 1321 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1322 } 1323 } 1324 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1325 PetscFunctionReturn(0); 1326 } 1327 1328 PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a) 1329 { 1330 PetscErrorCode ierr; 1331 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data; 1332 1333 PetscFunctionBegin; 1334 if (!Y->preallocated || !aij->nz) { 1335 ierr = MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);CHKERRQ(ierr); 1336 } 1337 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 1338 PetscFunctionReturn(0); 1339 } 1340 1341 /* -------------------------------------------------------------------*/ 1342 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1343 MatGetRow_SeqSBAIJ, 1344 MatRestoreRow_SeqSBAIJ, 1345 MatMult_SeqSBAIJ_N, 1346 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1347 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1348 MatMultAdd_SeqSBAIJ_N, 1349 0, 1350 0, 1351 0, 1352 /* 10*/ 0, 1353 0, 1354 MatCholeskyFactor_SeqSBAIJ, 1355 MatSOR_SeqSBAIJ, 1356 MatTranspose_SeqSBAIJ, 1357 /* 15*/ MatGetInfo_SeqSBAIJ, 1358 MatEqual_SeqSBAIJ, 1359 MatGetDiagonal_SeqSBAIJ, 1360 MatDiagonalScale_SeqSBAIJ, 1361 MatNorm_SeqSBAIJ, 1362 /* 20*/ 0, 1363 MatAssemblyEnd_SeqSBAIJ, 1364 MatSetOption_SeqSBAIJ, 1365 MatZeroEntries_SeqSBAIJ, 1366 /* 24*/ 0, 1367 0, 1368 0, 1369 0, 1370 0, 1371 /* 29*/ MatSetUp_SeqSBAIJ, 1372 0, 1373 0, 1374 0, 1375 0, 1376 /* 34*/ MatDuplicate_SeqSBAIJ, 1377 0, 1378 0, 1379 0, 1380 MatICCFactor_SeqSBAIJ, 1381 /* 39*/ MatAXPY_SeqSBAIJ, 1382 MatCreateSubMatrices_SeqSBAIJ, 1383 MatIncreaseOverlap_SeqSBAIJ, 1384 MatGetValues_SeqSBAIJ, 1385 MatCopy_SeqSBAIJ, 1386 /* 44*/ 0, 1387 MatScale_SeqSBAIJ, 1388 MatShift_SeqSBAIJ, 1389 0, 1390 MatZeroRowsColumns_SeqSBAIJ, 1391 /* 49*/ 0, 1392 MatGetRowIJ_SeqSBAIJ, 1393 MatRestoreRowIJ_SeqSBAIJ, 1394 0, 1395 0, 1396 /* 54*/ 0, 1397 0, 1398 0, 1399 0, 1400 MatSetValuesBlocked_SeqSBAIJ, 1401 /* 59*/ MatCreateSubMatrix_SeqSBAIJ, 1402 0, 1403 0, 1404 0, 1405 0, 1406 /* 64*/ 0, 1407 0, 1408 0, 1409 0, 1410 0, 1411 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1412 0, 1413 0, 1414 0, 1415 0, 1416 /* 74*/ 0, 1417 0, 1418 0, 1419 0, 1420 0, 1421 /* 79*/ 0, 1422 0, 1423 0, 1424 MatGetInertia_SeqSBAIJ, 1425 MatLoad_SeqSBAIJ, 1426 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1427 MatIsHermitian_SeqSBAIJ, 1428 MatIsStructurallySymmetric_SeqSBAIJ, 1429 0, 1430 0, 1431 /* 89*/ 0, 1432 0, 1433 0, 1434 0, 1435 0, 1436 /* 94*/ 0, 1437 0, 1438 0, 1439 0, 1440 0, 1441 /* 99*/ 0, 1442 0, 1443 0, 1444 0, 1445 0, 1446 /*104*/ 0, 1447 MatRealPart_SeqSBAIJ, 1448 MatImaginaryPart_SeqSBAIJ, 1449 MatGetRowUpperTriangular_SeqSBAIJ, 1450 MatRestoreRowUpperTriangular_SeqSBAIJ, 1451 /*109*/ 0, 1452 0, 1453 0, 1454 0, 1455 MatMissingDiagonal_SeqSBAIJ, 1456 /*114*/ 0, 1457 0, 1458 0, 1459 0, 1460 0, 1461 /*119*/ 0, 1462 0, 1463 0, 1464 0, 1465 0, 1466 /*124*/ 0, 1467 0, 1468 0, 1469 0, 1470 0, 1471 /*129*/ 0, 1472 0, 1473 0, 1474 0, 1475 0, 1476 /*134*/ 0, 1477 0, 1478 0, 1479 0, 1480 0, 1481 /*139*/ MatSetBlockSizes_Default, 1482 0, 1483 0, 1484 0, 1485 0, 1486 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ 1487 }; 1488 1489 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1490 { 1491 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1492 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1493 PetscErrorCode ierr; 1494 1495 PetscFunctionBegin; 1496 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1497 1498 /* allocate space for values if not already there */ 1499 if (!aij->saved_values) { 1500 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 1501 } 1502 1503 /* copy values over */ 1504 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1505 PetscFunctionReturn(0); 1506 } 1507 1508 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1509 { 1510 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1511 PetscErrorCode ierr; 1512 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1513 1514 PetscFunctionBegin; 1515 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1516 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1517 1518 /* copy values over */ 1519 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1520 PetscFunctionReturn(0); 1521 } 1522 1523 static PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1524 { 1525 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1526 PetscErrorCode ierr; 1527 PetscInt i,mbs,nbs,bs2; 1528 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1529 1530 PetscFunctionBegin; 1531 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1532 1533 ierr = MatSetBlockSize(B,PetscAbs(bs));CHKERRQ(ierr); 1534 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1535 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1536 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1537 1538 B->preallocated = PETSC_TRUE; 1539 1540 mbs = B->rmap->N/bs; 1541 nbs = B->cmap->n/bs; 1542 bs2 = bs*bs; 1543 1544 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"); 1545 1546 if (nz == MAT_SKIP_ALLOCATION) { 1547 skipallocation = PETSC_TRUE; 1548 nz = 0; 1549 } 1550 1551 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1552 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1553 if (nnz) { 1554 for (i=0; i<mbs; i++) { 1555 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]); 1556 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); 1557 } 1558 } 1559 1560 B->ops->mult = MatMult_SeqSBAIJ_N; 1561 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1562 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1563 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1564 1565 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1566 if (!flg) { 1567 switch (bs) { 1568 case 1: 1569 B->ops->mult = MatMult_SeqSBAIJ_1; 1570 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1571 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1572 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1573 break; 1574 case 2: 1575 B->ops->mult = MatMult_SeqSBAIJ_2; 1576 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1577 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1578 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1579 break; 1580 case 3: 1581 B->ops->mult = MatMult_SeqSBAIJ_3; 1582 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1583 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1584 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1585 break; 1586 case 4: 1587 B->ops->mult = MatMult_SeqSBAIJ_4; 1588 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1589 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1590 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1591 break; 1592 case 5: 1593 B->ops->mult = MatMult_SeqSBAIJ_5; 1594 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1595 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1596 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1597 break; 1598 case 6: 1599 B->ops->mult = MatMult_SeqSBAIJ_6; 1600 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1601 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1602 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1603 break; 1604 case 7: 1605 B->ops->mult = MatMult_SeqSBAIJ_7; 1606 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1607 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1608 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1609 break; 1610 } 1611 } 1612 1613 b->mbs = mbs; 1614 b->nbs = nbs; 1615 if (!skipallocation) { 1616 if (!b->imax) { 1617 ierr = PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);CHKERRQ(ierr); 1618 1619 b->free_imax_ilen = PETSC_TRUE; 1620 1621 ierr = PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1622 } 1623 if (!nnz) { 1624 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1625 else if (nz <= 0) nz = 1; 1626 for (i=0; i<mbs; i++) b->imax[i] = nz; 1627 nz = nz*mbs; /* total nz */ 1628 } else { 1629 nz = 0; 1630 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1631 } 1632 /* b->ilen will count nonzeros in each block row so far. */ 1633 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1634 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1635 1636 /* allocate the matrix space */ 1637 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1638 ierr = PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);CHKERRQ(ierr); 1639 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1640 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1641 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1642 1643 b->singlemalloc = PETSC_TRUE; 1644 1645 /* pointer to beginning of each row */ 1646 b->i[0] = 0; 1647 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1648 1649 b->free_a = PETSC_TRUE; 1650 b->free_ij = PETSC_TRUE; 1651 } else { 1652 b->free_a = PETSC_FALSE; 1653 b->free_ij = PETSC_FALSE; 1654 } 1655 1656 b->bs2 = bs2; 1657 b->nz = 0; 1658 b->maxnz = nz; 1659 b->inew = 0; 1660 b->jnew = 0; 1661 b->anew = 0; 1662 b->a2anew = 0; 1663 b->permute = PETSC_FALSE; 1664 1665 B->was_assembled = PETSC_FALSE; 1666 B->assembled = PETSC_FALSE; 1667 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1668 PetscFunctionReturn(0); 1669 } 1670 1671 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1672 { 1673 PetscInt i,j,m,nz,nz_max=0,*nnz; 1674 PetscScalar *values=0; 1675 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1676 PetscErrorCode ierr; 1677 PetscFunctionBegin; 1678 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs); 1679 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1680 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1681 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1682 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1683 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1684 m = B->rmap->n/bs; 1685 1686 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]); 1687 ierr = PetscMalloc1(m+1,&nnz);CHKERRQ(ierr); 1688 for (i=0; i<m; i++) { 1689 nz = ii[i+1] - ii[i]; 1690 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz); 1691 nz_max = PetscMax(nz_max,nz); 1692 nnz[i] = nz; 1693 } 1694 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1695 ierr = PetscFree(nnz);CHKERRQ(ierr); 1696 1697 values = (PetscScalar*)V; 1698 if (!values) { 1699 ierr = PetscCalloc1(bs*bs*nz_max,&values);CHKERRQ(ierr); 1700 } 1701 for (i=0; i<m; i++) { 1702 PetscInt ncols = ii[i+1] - ii[i]; 1703 const PetscInt *icols = jj + ii[i]; 1704 if (!roworiented || bs == 1) { 1705 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1706 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr); 1707 } else { 1708 for (j=0; j<ncols; j++) { 1709 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1710 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr); 1711 } 1712 } 1713 } 1714 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); } 1715 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1716 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1717 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1718 PetscFunctionReturn(0); 1719 } 1720 1721 /* 1722 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1723 */ 1724 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1725 { 1726 PetscErrorCode ierr; 1727 PetscBool flg = PETSC_FALSE; 1728 PetscInt bs = B->rmap->bs; 1729 1730 PetscFunctionBegin; 1731 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1732 if (flg) bs = 8; 1733 1734 if (!natural) { 1735 switch (bs) { 1736 case 1: 1737 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1738 break; 1739 case 2: 1740 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1741 break; 1742 case 3: 1743 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1744 break; 1745 case 4: 1746 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1747 break; 1748 case 5: 1749 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1750 break; 1751 case 6: 1752 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1753 break; 1754 case 7: 1755 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1756 break; 1757 default: 1758 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1759 break; 1760 } 1761 } else { 1762 switch (bs) { 1763 case 1: 1764 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1765 break; 1766 case 2: 1767 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1768 break; 1769 case 3: 1770 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1771 break; 1772 case 4: 1773 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1774 break; 1775 case 5: 1776 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1777 break; 1778 case 6: 1779 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1780 break; 1781 case 7: 1782 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1783 break; 1784 default: 1785 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1786 break; 1787 } 1788 } 1789 PetscFunctionReturn(0); 1790 } 1791 1792 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1793 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1794 1795 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1796 { 1797 PetscInt n = A->rmap->n; 1798 PetscErrorCode ierr; 1799 1800 PetscFunctionBegin; 1801 #if defined(PETSC_USE_COMPLEX) 1802 if (A->hermitian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported"); 1803 #endif 1804 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 1805 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1806 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1807 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1808 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1809 1810 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1811 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1812 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1813 1814 (*B)->factortype = ftype; 1815 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 1816 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 1817 PetscFunctionReturn(0); 1818 } 1819 1820 /*@C 1821 MatSeqSBAIJGetArray - gives access to the array where the data for a MATSEQSBAIJ matrix is stored 1822 1823 Not Collective 1824 1825 Input Parameter: 1826 . mat - a MATSEQSBAIJ matrix 1827 1828 Output Parameter: 1829 . array - pointer to the data 1830 1831 Level: intermediate 1832 1833 .seealso: MatSeqSBAIJRestoreArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray() 1834 @*/ 1835 PetscErrorCode MatSeqSBAIJGetArray(Mat A,PetscScalar **array) 1836 { 1837 PetscErrorCode ierr; 1838 1839 PetscFunctionBegin; 1840 ierr = PetscUseMethod(A,"MatSeqSBAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1841 PetscFunctionReturn(0); 1842 } 1843 1844 /*@C 1845 MatSeqSBAIJRestoreArray - returns access to the array where the data for a MATSEQSBAIJ matrix is stored obtained by MatSeqSBAIJGetArray() 1846 1847 Not Collective 1848 1849 Input Parameters: 1850 . mat - a MATSEQSBAIJ matrix 1851 . array - pointer to the data 1852 1853 Level: intermediate 1854 1855 .seealso: MatSeqSBAIJGetArray(), MatSeqAIJGetArray(), MatSeqAIJRestoreArray() 1856 @*/ 1857 PetscErrorCode MatSeqSBAIJRestoreArray(Mat A,PetscScalar **array) 1858 { 1859 PetscErrorCode ierr; 1860 1861 PetscFunctionBegin; 1862 ierr = PetscUseMethod(A,"MatSeqSBAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1863 PetscFunctionReturn(0); 1864 } 1865 1866 /*MC 1867 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1868 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1869 1870 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1871 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1872 1873 Options Database Keys: 1874 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1875 1876 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1877 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1878 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. 1879 1880 1881 Level: beginner 1882 1883 .seealso: MatCreateSeqSBAIJ 1884 M*/ 1885 1886 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1887 1888 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1889 { 1890 Mat_SeqSBAIJ *b; 1891 PetscErrorCode ierr; 1892 PetscMPIInt size; 1893 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1894 1895 PetscFunctionBegin; 1896 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 1897 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1898 1899 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1900 B->data = (void*)b; 1901 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1902 1903 B->ops->destroy = MatDestroy_SeqSBAIJ; 1904 B->ops->view = MatView_SeqSBAIJ; 1905 b->row = 0; 1906 b->icol = 0; 1907 b->reallocs = 0; 1908 b->saved_values = 0; 1909 b->inode.limit = 5; 1910 b->inode.max_limit = 5; 1911 1912 b->roworiented = PETSC_TRUE; 1913 b->nonew = 0; 1914 b->diag = 0; 1915 b->solve_work = 0; 1916 b->mult_work = 0; 1917 B->spptr = 0; 1918 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1919 b->keepnonzeropattern = PETSC_FALSE; 1920 1921 b->inew = 0; 1922 b->jnew = 0; 1923 b->anew = 0; 1924 b->a2anew = 0; 1925 b->permute = PETSC_FALSE; 1926 1927 b->ignore_ltriangular = PETSC_TRUE; 1928 1929 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1930 1931 b->getrow_utriangular = PETSC_FALSE; 1932 1933 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 1934 1935 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJGetArray_C",MatSeqSBAIJGetArray_SeqSBAIJ);CHKERRQ(ierr); 1936 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJRestoreArray_C",MatSeqSBAIJRestoreArray_SeqSBAIJ);CHKERRQ(ierr); 1937 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1938 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1939 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1940 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1941 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1942 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1943 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr); 1944 #if defined(PETSC_HAVE_ELEMENTAL) 1945 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);CHKERRQ(ierr); 1946 #endif 1947 1948 B->symmetric = PETSC_TRUE; 1949 B->structurally_symmetric = PETSC_TRUE; 1950 B->symmetric_set = PETSC_TRUE; 1951 B->structurally_symmetric_set = PETSC_TRUE; 1952 1953 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1954 1955 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1956 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 1957 if (no_unroll) { 1958 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 1959 } 1960 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 1961 if (no_inode) { 1962 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 1963 } 1964 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 1965 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1966 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1967 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1968 PetscFunctionReturn(0); 1969 } 1970 1971 /*@C 1972 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1973 compressed row) format. For good matrix assembly performance the 1974 user should preallocate the matrix storage by setting the parameter nz 1975 (or the array nnz). By setting these parameters accurately, performance 1976 during matrix assembly can be increased by more than a factor of 50. 1977 1978 Collective on Mat 1979 1980 Input Parameters: 1981 + B - the symmetric matrix 1982 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 1983 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 1984 . nz - number of block nonzeros per block row (same for all rows) 1985 - nnz - array containing the number of block nonzeros in the upper triangular plus 1986 diagonal portion of each block (possibly different for each block row) or NULL 1987 1988 Options Database Keys: 1989 . -mat_no_unroll - uses code that does not unroll the loops in the 1990 block calculations (much slower) 1991 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1992 1993 Level: intermediate 1994 1995 Notes: 1996 Specify the preallocated storage with either nz or nnz (not both). 1997 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 1998 allocation. See Users-Manual: ch_mat for details. 1999 2000 You can call MatGetInfo() to get information on how effective the preallocation was; 2001 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2002 You can also run with the option -info and look for messages with the string 2003 malloc in them to see if additional memory allocation was needed. 2004 2005 If the nnz parameter is given then the nz parameter is ignored 2006 2007 2008 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2009 @*/ 2010 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2011 { 2012 PetscErrorCode ierr; 2013 2014 PetscFunctionBegin; 2015 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2016 PetscValidType(B,1); 2017 PetscValidLogicalCollectiveInt(B,bs,2); 2018 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2019 PetscFunctionReturn(0); 2020 } 2021 2022 /*@C 2023 MatSeqSBAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in symmetric block AIJ format. 2024 2025 Input Parameters: 2026 + B - the matrix 2027 . bs - size of block, the blocks are ALWAYS square. 2028 . i - the indices into j for the start of each local row (starts with zero) 2029 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2030 - v - optional values in the matrix 2031 2032 Level: developer 2033 2034 Notes: 2035 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2036 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2037 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2038 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2039 block column and the second index is over columns within a block. 2040 2041 .keywords: matrix, block, aij, compressed row, sparse 2042 2043 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ 2044 @*/ 2045 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2046 { 2047 PetscErrorCode ierr; 2048 2049 PetscFunctionBegin; 2050 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2051 PetscValidType(B,1); 2052 PetscValidLogicalCollectiveInt(B,bs,2); 2053 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2054 PetscFunctionReturn(0); 2055 } 2056 2057 /*@C 2058 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2059 compressed row) format. For good matrix assembly performance the 2060 user should preallocate the matrix storage by setting the parameter nz 2061 (or the array nnz). By setting these parameters accurately, performance 2062 during matrix assembly can be increased by more than a factor of 50. 2063 2064 Collective on MPI_Comm 2065 2066 Input Parameters: 2067 + comm - MPI communicator, set to PETSC_COMM_SELF 2068 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2069 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2070 . m - number of rows, or number of columns 2071 . nz - number of block nonzeros per block row (same for all rows) 2072 - nnz - array containing the number of block nonzeros in the upper triangular plus 2073 diagonal portion of each block (possibly different for each block row) or NULL 2074 2075 Output Parameter: 2076 . A - the symmetric matrix 2077 2078 Options Database Keys: 2079 . -mat_no_unroll - uses code that does not unroll the loops in the 2080 block calculations (much slower) 2081 . -mat_block_size - size of the blocks to use 2082 2083 Level: intermediate 2084 2085 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2086 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2087 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2088 2089 Notes: 2090 The number of rows and columns must be divisible by blocksize. 2091 This matrix type does not support complex Hermitian operation. 2092 2093 Specify the preallocated storage with either nz or nnz (not both). 2094 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2095 allocation. See Users-Manual: ch_mat for details. 2096 2097 If the nnz parameter is given then the nz parameter is ignored 2098 2099 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2100 @*/ 2101 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2102 { 2103 PetscErrorCode ierr; 2104 2105 PetscFunctionBegin; 2106 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2107 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2108 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2109 ierr = MatSeqSBAIJSetPreallocation(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2110 PetscFunctionReturn(0); 2111 } 2112 2113 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2114 { 2115 Mat C; 2116 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2117 PetscErrorCode ierr; 2118 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2119 2120 PetscFunctionBegin; 2121 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2122 2123 *B = 0; 2124 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2125 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2126 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2127 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2128 c = (Mat_SeqSBAIJ*)C->data; 2129 2130 C->preallocated = PETSC_TRUE; 2131 C->factortype = A->factortype; 2132 c->row = 0; 2133 c->icol = 0; 2134 c->saved_values = 0; 2135 c->keepnonzeropattern = a->keepnonzeropattern; 2136 C->assembled = PETSC_TRUE; 2137 2138 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2139 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2140 c->bs2 = a->bs2; 2141 c->mbs = a->mbs; 2142 c->nbs = a->nbs; 2143 2144 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2145 c->imax = a->imax; 2146 c->ilen = a->ilen; 2147 c->free_imax_ilen = PETSC_FALSE; 2148 } else { 2149 ierr = PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);CHKERRQ(ierr); 2150 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2151 for (i=0; i<mbs; i++) { 2152 c->imax[i] = a->imax[i]; 2153 c->ilen[i] = a->ilen[i]; 2154 } 2155 c->free_imax_ilen = PETSC_TRUE; 2156 } 2157 2158 /* allocate the matrix space */ 2159 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2160 ierr = PetscMalloc1(bs2*nz,&c->a);CHKERRQ(ierr); 2161 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2162 c->i = a->i; 2163 c->j = a->j; 2164 c->singlemalloc = PETSC_FALSE; 2165 c->free_a = PETSC_TRUE; 2166 c->free_ij = PETSC_FALSE; 2167 c->parent = A; 2168 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2169 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2170 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2171 } else { 2172 ierr = PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);CHKERRQ(ierr); 2173 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2174 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2175 c->singlemalloc = PETSC_TRUE; 2176 c->free_a = PETSC_TRUE; 2177 c->free_ij = PETSC_TRUE; 2178 } 2179 if (mbs > 0) { 2180 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2181 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2182 } 2183 if (cpvalues == MAT_COPY_VALUES) { 2184 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2185 } else { 2186 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2187 } 2188 if (a->jshort) { 2189 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2190 /* if the parent matrix is reassembled, this child matrix will never notice */ 2191 ierr = PetscMalloc1(nz,&c->jshort);CHKERRQ(ierr); 2192 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2193 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2194 2195 c->free_jshort = PETSC_TRUE; 2196 } 2197 } 2198 2199 c->roworiented = a->roworiented; 2200 c->nonew = a->nonew; 2201 2202 if (a->diag) { 2203 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2204 c->diag = a->diag; 2205 c->free_diag = PETSC_FALSE; 2206 } else { 2207 ierr = PetscMalloc1(mbs,&c->diag);CHKERRQ(ierr); 2208 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2209 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2210 c->free_diag = PETSC_TRUE; 2211 } 2212 } 2213 c->nz = a->nz; 2214 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2215 c->solve_work = 0; 2216 c->mult_work = 0; 2217 2218 *B = C; 2219 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2220 PetscFunctionReturn(0); 2221 } 2222 2223 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2224 { 2225 Mat_SeqSBAIJ *a; 2226 PetscErrorCode ierr; 2227 int fd; 2228 PetscMPIInt size; 2229 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs = newmat->rmap->bs; 2230 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2231 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2232 PetscInt *masked,nmask,tmp,bs2,ishift; 2233 PetscScalar *aa; 2234 MPI_Comm comm; 2235 2236 PetscFunctionBegin; 2237 /* force binary viewer to load .info file if it has not yet done so */ 2238 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 2239 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2240 ierr = PetscOptionsGetInt(((PetscObject)newmat)->options,((PetscObject)newmat)->prefix,"-matload_block_size",&bs,NULL);CHKERRQ(ierr); 2241 if (bs < 0) bs = 1; 2242 bs2 = bs*bs; 2243 2244 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2245 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2246 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2247 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2248 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2249 M = header[1]; N = header[2]; nz = header[3]; 2250 2251 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2252 2253 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2254 2255 /* 2256 This code adds extra rows to make sure the number of rows is 2257 divisible by the blocksize 2258 */ 2259 mbs = M/bs; 2260 extra_rows = bs - M + bs*(mbs); 2261 if (extra_rows == bs) extra_rows = 0; 2262 else mbs++; 2263 if (extra_rows) { 2264 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2265 } 2266 2267 /* Set global sizes if not already set */ 2268 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2269 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2270 } else { /* Check if the matrix global sizes are correct */ 2271 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2272 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); 2273 } 2274 2275 /* read in row lengths */ 2276 ierr = PetscMalloc1(M+extra_rows,&rowlengths);CHKERRQ(ierr); 2277 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2278 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2279 2280 /* read in column indices */ 2281 ierr = PetscMalloc1(nz+extra_rows,&jj);CHKERRQ(ierr); 2282 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2283 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2284 2285 /* loop over row lengths determining block row lengths */ 2286 ierr = PetscCalloc1(mbs,&s_browlengths);CHKERRQ(ierr); 2287 ierr = PetscMalloc2(mbs,&mask,mbs,&masked);CHKERRQ(ierr); 2288 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2289 rowcount = 0; 2290 nzcount = 0; 2291 for (i=0; i<mbs; i++) { 2292 nmask = 0; 2293 for (j=0; j<bs; j++) { 2294 kmax = rowlengths[rowcount]; 2295 for (k=0; k<kmax; k++) { 2296 tmp = jj[nzcount++]/bs; /* block col. index */ 2297 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2298 } 2299 rowcount++; 2300 } 2301 s_browlengths[i] += nmask; 2302 2303 /* zero out the mask elements we set */ 2304 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2305 } 2306 2307 /* Do preallocation */ 2308 ierr = MatSeqSBAIJSetPreallocation(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2309 a = (Mat_SeqSBAIJ*)newmat->data; 2310 2311 /* set matrix "i" values */ 2312 a->i[0] = 0; 2313 for (i=1; i<= mbs; i++) { 2314 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2315 a->ilen[i-1] = s_browlengths[i-1]; 2316 } 2317 a->nz = a->i[mbs]; 2318 2319 /* read in nonzero values */ 2320 ierr = PetscMalloc1(nz+extra_rows,&aa);CHKERRQ(ierr); 2321 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2322 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2323 2324 /* set "a" and "j" values into matrix */ 2325 nzcount = 0; jcount = 0; 2326 for (i=0; i<mbs; i++) { 2327 nzcountb = nzcount; 2328 nmask = 0; 2329 for (j=0; j<bs; j++) { 2330 kmax = rowlengths[i*bs+j]; 2331 for (k=0; k<kmax; k++) { 2332 tmp = jj[nzcount++]/bs; /* block col. index */ 2333 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2334 } 2335 } 2336 /* sort the masked values */ 2337 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2338 2339 /* set "j" values into matrix */ 2340 maskcount = 1; 2341 for (j=0; j<nmask; j++) { 2342 a->j[jcount++] = masked[j]; 2343 mask[masked[j]] = maskcount++; 2344 } 2345 2346 /* set "a" values into matrix */ 2347 ishift = bs2*a->i[i]; 2348 for (j=0; j<bs; j++) { 2349 kmax = rowlengths[i*bs+j]; 2350 for (k=0; k<kmax; k++) { 2351 tmp = jj[nzcountb]/bs; /* block col. index */ 2352 if (tmp >= i) { 2353 block = mask[tmp] - 1; 2354 point = jj[nzcountb] - bs*tmp; 2355 idx = ishift + bs2*block + j + bs*point; 2356 a->a[idx] = aa[nzcountb]; 2357 } 2358 nzcountb++; 2359 } 2360 } 2361 /* zero out the mask elements we set */ 2362 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2363 } 2364 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2365 2366 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2367 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2368 ierr = PetscFree(aa);CHKERRQ(ierr); 2369 ierr = PetscFree(jj);CHKERRQ(ierr); 2370 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2371 2372 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2373 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2374 PetscFunctionReturn(0); 2375 } 2376 2377 /*@ 2378 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2379 (upper triangular entries in CSR format) provided by the user. 2380 2381 Collective on MPI_Comm 2382 2383 Input Parameters: 2384 + comm - must be an MPI communicator of size 1 2385 . bs - size of block 2386 . m - number of rows 2387 . n - number of columns 2388 . i - row indices 2389 . j - column indices 2390 - a - matrix values 2391 2392 Output Parameter: 2393 . mat - the matrix 2394 2395 Level: advanced 2396 2397 Notes: 2398 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2399 once the matrix is destroyed 2400 2401 You cannot set new nonzero locations into this matrix, that will generate an error. 2402 2403 The i and j indices are 0 based 2404 2405 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 2406 it is the regular CSR format excluding the lower triangular elements. 2407 2408 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2409 2410 @*/ 2411 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 2412 { 2413 PetscErrorCode ierr; 2414 PetscInt ii; 2415 Mat_SeqSBAIJ *sbaij; 2416 2417 PetscFunctionBegin; 2418 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2419 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2420 2421 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2422 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2423 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2424 ierr = MatSeqSBAIJSetPreallocation(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2425 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2426 ierr = PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);CHKERRQ(ierr); 2427 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2428 2429 sbaij->i = i; 2430 sbaij->j = j; 2431 sbaij->a = a; 2432 2433 sbaij->singlemalloc = PETSC_FALSE; 2434 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2435 sbaij->free_a = PETSC_FALSE; 2436 sbaij->free_ij = PETSC_FALSE; 2437 2438 for (ii=0; ii<m; ii++) { 2439 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2440 #if defined(PETSC_USE_DEBUG) 2441 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]); 2442 #endif 2443 } 2444 #if defined(PETSC_USE_DEBUG) 2445 for (ii=0; ii<sbaij->i[m]; ii++) { 2446 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2447 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]); 2448 } 2449 #endif 2450 2451 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2452 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2453 PetscFunctionReturn(0); 2454 } 2455 2456 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2457 { 2458 PetscErrorCode ierr; 2459 PetscMPIInt size; 2460 2461 PetscFunctionBegin; 2462 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2463 if (size == 1 && scall == MAT_REUSE_MATRIX) { 2464 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2465 } else { 2466 ierr = MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2467 } 2468 PetscFunctionReturn(0); 2469 } 2470 2471 2472 2473 2474