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