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