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