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