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