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