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