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 1023 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1024 ierr = MatSeqSBAIJSetNumericFactorization_inplace(inA,row_identity);CHKERRQ(ierr); 1025 1026 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1027 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1028 a->row = row; 1029 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1030 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1031 a->col = row; 1032 1033 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1034 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1035 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 1036 1037 if (!a->solve_work) { 1038 ierr = PetscMalloc1(inA->rmap->N+inA->rmap->bs,&a->solve_work);CHKERRQ(ierr); 1039 ierr = PetscLogObjectMemory((PetscObject)inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1040 } 1041 1042 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1043 PetscFunctionReturn(0); 1044 } 1045 1046 #undef __FUNCT__ 1047 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1048 PetscErrorCode MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1049 { 1050 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ*)mat->data; 1051 PetscInt i,nz,n; 1052 PetscErrorCode ierr; 1053 1054 PetscFunctionBegin; 1055 nz = baij->maxnz; 1056 n = mat->cmap->n; 1057 for (i=0; i<nz; i++) baij->j[i] = indices[i]; 1058 1059 baij->nz = nz; 1060 for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i]; 1061 1062 ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1063 PetscFunctionReturn(0); 1064 } 1065 1066 #undef __FUNCT__ 1067 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1068 /*@ 1069 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1070 in the matrix. 1071 1072 Input Parameters: 1073 + mat - the SeqSBAIJ matrix 1074 - indices - the column indices 1075 1076 Level: advanced 1077 1078 Notes: 1079 This can be called if you have precomputed the nonzero structure of the 1080 matrix and want to provide it to the matrix object to improve the performance 1081 of the MatSetValues() operation. 1082 1083 You MUST have set the correct numbers of nonzeros per row in the call to 1084 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1085 1086 MUST be called before any calls to MatSetValues() 1087 1088 .seealso: MatCreateSeqSBAIJ 1089 @*/ 1090 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1091 { 1092 PetscErrorCode ierr; 1093 1094 PetscFunctionBegin; 1095 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1096 PetscValidPointer(indices,2); 1097 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 1098 PetscFunctionReturn(0); 1099 } 1100 1101 #undef __FUNCT__ 1102 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1103 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1104 { 1105 PetscErrorCode ierr; 1106 1107 PetscFunctionBegin; 1108 /* If the two matrices have the same copy implementation, use fast copy. */ 1109 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1110 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1111 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1112 1113 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"); 1114 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1115 } else { 1116 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1117 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1118 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1119 } 1120 PetscFunctionReturn(0); 1121 } 1122 1123 #undef __FUNCT__ 1124 #define __FUNCT__ "MatSetUp_SeqSBAIJ" 1125 PetscErrorCode MatSetUp_SeqSBAIJ(Mat A) 1126 { 1127 PetscErrorCode ierr; 1128 1129 PetscFunctionBegin; 1130 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,A->rmap->bs,PETSC_DEFAULT,0);CHKERRQ(ierr); 1131 PetscFunctionReturn(0); 1132 } 1133 1134 #undef __FUNCT__ 1135 #define __FUNCT__ "MatSeqSBAIJGetArray_SeqSBAIJ" 1136 PetscErrorCode MatSeqSBAIJGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1137 { 1138 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1139 1140 PetscFunctionBegin; 1141 *array = a->a; 1142 PetscFunctionReturn(0); 1143 } 1144 1145 #undef __FUNCT__ 1146 #define __FUNCT__ "MatSeqSBAIJRestoreArray_SeqSBAIJ" 1147 PetscErrorCode MatSeqSBAIJRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1148 { 1149 PetscFunctionBegin; 1150 PetscFunctionReturn(0); 1151 } 1152 1153 #undef __FUNCT__ 1154 #define __FUNCT__ "MatAXPYGetPreallocation_SeqSBAIJ" 1155 PetscErrorCode MatAXPYGetPreallocation_SeqSBAIJ(Mat Y,Mat X,PetscInt *nnz) 1156 { 1157 PetscInt bs = Y->rmap->bs,mbs = Y->rmap->N/bs; 1158 Mat_SeqSBAIJ *x = (Mat_SeqSBAIJ*)X->data; 1159 Mat_SeqSBAIJ *y = (Mat_SeqSBAIJ*)Y->data; 1160 PetscErrorCode ierr; 1161 1162 PetscFunctionBegin; 1163 /* Set the number of nonzeros in the new matrix */ 1164 ierr = MatAXPYGetPreallocation_SeqX_private(mbs,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 1165 PetscFunctionReturn(0); 1166 } 1167 1168 #undef __FUNCT__ 1169 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1170 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1171 { 1172 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ*)X->data, *y=(Mat_SeqSBAIJ*)Y->data; 1173 PetscErrorCode ierr; 1174 PetscInt bs=Y->rmap->bs,bs2=bs*bs; 1175 PetscBLASInt one = 1; 1176 1177 PetscFunctionBegin; 1178 if (str == SAME_NONZERO_PATTERN) { 1179 PetscScalar alpha = a; 1180 PetscBLASInt bnz; 1181 ierr = PetscBLASIntCast(x->nz*bs2,&bnz);CHKERRQ(ierr); 1182 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 1183 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 1184 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1185 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE);CHKERRQ(ierr); 1186 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1187 ierr = MatSetOption(X,MAT_GETROW_UPPERTRIANGULAR,PETSC_FALSE);CHKERRQ(ierr); 1188 } else { 1189 Mat B; 1190 PetscInt *nnz; 1191 if (bs != X->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrices must have same block size"); 1192 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1193 ierr = MatGetRowUpperTriangular(Y);CHKERRQ(ierr); 1194 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 1195 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 1196 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 1197 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 1198 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 1199 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 1200 ierr = MatAXPYGetPreallocation_SeqSBAIJ(Y,X,nnz);CHKERRQ(ierr); 1201 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1202 1203 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 1204 1205 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 1206 ierr = PetscFree(nnz);CHKERRQ(ierr); 1207 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1208 ierr = MatRestoreRowUpperTriangular(Y);CHKERRQ(ierr); 1209 } 1210 PetscFunctionReturn(0); 1211 } 1212 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1215 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1216 { 1217 PetscFunctionBegin; 1218 *flg = PETSC_TRUE; 1219 PetscFunctionReturn(0); 1220 } 1221 1222 #undef __FUNCT__ 1223 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1224 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1225 { 1226 PetscFunctionBegin; 1227 *flg = PETSC_TRUE; 1228 PetscFunctionReturn(0); 1229 } 1230 1231 #undef __FUNCT__ 1232 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1233 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1234 { 1235 PetscFunctionBegin; 1236 *flg = PETSC_FALSE; 1237 PetscFunctionReturn(0); 1238 } 1239 1240 #undef __FUNCT__ 1241 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1242 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1243 { 1244 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1245 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1246 MatScalar *aa = a->a; 1247 1248 PetscFunctionBegin; 1249 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1250 PetscFunctionReturn(0); 1251 } 1252 1253 #undef __FUNCT__ 1254 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1255 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1256 { 1257 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1258 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1259 MatScalar *aa = a->a; 1260 1261 PetscFunctionBegin; 1262 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1263 PetscFunctionReturn(0); 1264 } 1265 1266 #undef __FUNCT__ 1267 #define __FUNCT__ "MatZeroRowsColumns_SeqSBAIJ" 1268 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1269 { 1270 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1271 PetscErrorCode ierr; 1272 PetscInt i,j,k,count; 1273 PetscInt bs =A->rmap->bs,bs2=baij->bs2,row,col; 1274 PetscScalar zero = 0.0; 1275 MatScalar *aa; 1276 const PetscScalar *xx; 1277 PetscScalar *bb; 1278 PetscBool *zeroed,vecs = PETSC_FALSE; 1279 1280 PetscFunctionBegin; 1281 /* fix right hand side if needed */ 1282 if (x && b) { 1283 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1284 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1285 vecs = PETSC_TRUE; 1286 } 1287 1288 /* zero the columns */ 1289 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1290 for (i=0; i<is_n; i++) { 1291 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]); 1292 zeroed[is_idx[i]] = PETSC_TRUE; 1293 } 1294 if (vecs) { 1295 for (i=0; i<A->rmap->N; i++) { 1296 row = i/bs; 1297 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1298 for (k=0; k<bs; k++) { 1299 col = bs*baij->j[j] + k; 1300 if (col <= i) continue; 1301 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1302 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col]; 1303 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i]; 1304 } 1305 } 1306 } 1307 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]]; 1308 } 1309 1310 for (i=0; i<A->rmap->N; i++) { 1311 if (!zeroed[i]) { 1312 row = i/bs; 1313 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1314 for (k=0; k<bs; k++) { 1315 col = bs*baij->j[j] + k; 1316 if (zeroed[col]) { 1317 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1318 aa[0] = 0.0; 1319 } 1320 } 1321 } 1322 } 1323 } 1324 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1325 if (vecs) { 1326 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1327 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1328 } 1329 1330 /* zero the rows */ 1331 for (i=0; i<is_n; i++) { 1332 row = is_idx[i]; 1333 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1334 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1335 for (k=0; k<count; k++) { 1336 aa[0] = zero; 1337 aa += bs; 1338 } 1339 if (diag != 0.0) { 1340 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1341 } 1342 } 1343 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1344 PetscFunctionReturn(0); 1345 } 1346 1347 #undef __FUNCT__ 1348 #define __FUNCT__ "MatShift_SeqSBAIJ" 1349 PetscErrorCode MatShift_SeqSBAIJ(Mat Y,PetscScalar a) 1350 { 1351 PetscErrorCode ierr; 1352 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)Y->data; 1353 1354 PetscFunctionBegin; 1355 if (!Y->preallocated || !aij->nz) { 1356 ierr = MatSeqSBAIJSetPreallocation(Y,Y->rmap->bs,1,NULL);CHKERRQ(ierr); 1357 } 1358 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 1359 PetscFunctionReturn(0); 1360 } 1361 1362 /* -------------------------------------------------------------------*/ 1363 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1364 MatGetRow_SeqSBAIJ, 1365 MatRestoreRow_SeqSBAIJ, 1366 MatMult_SeqSBAIJ_N, 1367 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1368 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1369 MatMultAdd_SeqSBAIJ_N, 1370 0, 1371 0, 1372 0, 1373 /* 10*/ 0, 1374 0, 1375 MatCholeskyFactor_SeqSBAIJ, 1376 MatSOR_SeqSBAIJ, 1377 MatTranspose_SeqSBAIJ, 1378 /* 15*/ MatGetInfo_SeqSBAIJ, 1379 MatEqual_SeqSBAIJ, 1380 MatGetDiagonal_SeqSBAIJ, 1381 MatDiagonalScale_SeqSBAIJ, 1382 MatNorm_SeqSBAIJ, 1383 /* 20*/ 0, 1384 MatAssemblyEnd_SeqSBAIJ, 1385 MatSetOption_SeqSBAIJ, 1386 MatZeroEntries_SeqSBAIJ, 1387 /* 24*/ 0, 1388 0, 1389 0, 1390 0, 1391 0, 1392 /* 29*/ MatSetUp_SeqSBAIJ, 1393 0, 1394 0, 1395 0, 1396 0, 1397 /* 34*/ MatDuplicate_SeqSBAIJ, 1398 0, 1399 0, 1400 0, 1401 MatICCFactor_SeqSBAIJ, 1402 /* 39*/ MatAXPY_SeqSBAIJ, 1403 MatGetSubMatrices_SeqSBAIJ, 1404 MatIncreaseOverlap_SeqSBAIJ, 1405 MatGetValues_SeqSBAIJ, 1406 MatCopy_SeqSBAIJ, 1407 /* 44*/ 0, 1408 MatScale_SeqSBAIJ, 1409 MatShift_SeqSBAIJ, 1410 0, 1411 MatZeroRowsColumns_SeqSBAIJ, 1412 /* 49*/ 0, 1413 MatGetRowIJ_SeqSBAIJ, 1414 MatRestoreRowIJ_SeqSBAIJ, 1415 0, 1416 0, 1417 /* 54*/ 0, 1418 0, 1419 0, 1420 0, 1421 MatSetValuesBlocked_SeqSBAIJ, 1422 /* 59*/ MatGetSubMatrix_SeqSBAIJ, 1423 0, 1424 0, 1425 0, 1426 0, 1427 /* 64*/ 0, 1428 0, 1429 0, 1430 0, 1431 0, 1432 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1433 0, 1434 0, 1435 0, 1436 0, 1437 /* 74*/ 0, 1438 0, 1439 0, 1440 0, 1441 0, 1442 /* 79*/ 0, 1443 0, 1444 0, 1445 MatGetInertia_SeqSBAIJ, 1446 MatLoad_SeqSBAIJ, 1447 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1448 MatIsHermitian_SeqSBAIJ, 1449 MatIsStructurallySymmetric_SeqSBAIJ, 1450 0, 1451 0, 1452 /* 89*/ 0, 1453 0, 1454 0, 1455 0, 1456 0, 1457 /* 94*/ 0, 1458 0, 1459 0, 1460 0, 1461 0, 1462 /* 99*/ 0, 1463 0, 1464 0, 1465 0, 1466 0, 1467 /*104*/ 0, 1468 MatRealPart_SeqSBAIJ, 1469 MatImaginaryPart_SeqSBAIJ, 1470 MatGetRowUpperTriangular_SeqSBAIJ, 1471 MatRestoreRowUpperTriangular_SeqSBAIJ, 1472 /*109*/ 0, 1473 0, 1474 0, 1475 0, 1476 MatMissingDiagonal_SeqSBAIJ, 1477 /*114*/ 0, 1478 0, 1479 0, 1480 0, 1481 0, 1482 /*119*/ 0, 1483 0, 1484 0, 1485 0, 1486 0, 1487 /*124*/ 0, 1488 0, 1489 0, 1490 0, 1491 0, 1492 /*129*/ 0, 1493 0, 1494 0, 1495 0, 1496 0, 1497 /*134*/ 0, 1498 0, 1499 0, 1500 0, 1501 0, 1502 /*139*/ 0, 1503 0, 1504 0, 1505 0, 1506 0, 1507 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ 1508 }; 1509 1510 #undef __FUNCT__ 1511 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1512 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1513 { 1514 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1515 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1516 PetscErrorCode ierr; 1517 1518 PetscFunctionBegin; 1519 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1520 1521 /* allocate space for values if not already there */ 1522 if (!aij->saved_values) { 1523 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 1524 } 1525 1526 /* copy values over */ 1527 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1528 PetscFunctionReturn(0); 1529 } 1530 1531 #undef __FUNCT__ 1532 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1533 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1534 { 1535 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1536 PetscErrorCode ierr; 1537 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1538 1539 PetscFunctionBegin; 1540 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1541 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1542 1543 /* copy values over */ 1544 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1545 PetscFunctionReturn(0); 1546 } 1547 1548 #undef __FUNCT__ 1549 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1550 PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1551 { 1552 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1553 PetscErrorCode ierr; 1554 PetscInt i,mbs,nbs,bs2; 1555 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1556 1557 PetscFunctionBegin; 1558 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1559 B->preallocated = PETSC_TRUE; 1560 1561 ierr = MatSetBlockSize(B,PetscAbs(bs));CHKERRQ(ierr); 1562 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1563 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1564 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1565 1566 mbs = B->rmap->N/bs; 1567 nbs = B->cmap->n/bs; 1568 bs2 = bs*bs; 1569 1570 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"); 1571 1572 if (nz == MAT_SKIP_ALLOCATION) { 1573 skipallocation = PETSC_TRUE; 1574 nz = 0; 1575 } 1576 1577 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1578 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1579 if (nnz) { 1580 for (i=0; i<mbs; i++) { 1581 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]); 1582 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); 1583 } 1584 } 1585 1586 B->ops->mult = MatMult_SeqSBAIJ_N; 1587 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1588 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1589 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1590 1591 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1592 if (!flg) { 1593 switch (bs) { 1594 case 1: 1595 B->ops->mult = MatMult_SeqSBAIJ_1; 1596 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1597 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1598 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1599 break; 1600 case 2: 1601 B->ops->mult = MatMult_SeqSBAIJ_2; 1602 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1603 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1604 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1605 break; 1606 case 3: 1607 B->ops->mult = MatMult_SeqSBAIJ_3; 1608 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1609 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1610 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1611 break; 1612 case 4: 1613 B->ops->mult = MatMult_SeqSBAIJ_4; 1614 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1615 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1616 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1617 break; 1618 case 5: 1619 B->ops->mult = MatMult_SeqSBAIJ_5; 1620 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1621 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1622 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1623 break; 1624 case 6: 1625 B->ops->mult = MatMult_SeqSBAIJ_6; 1626 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1627 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1628 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1629 break; 1630 case 7: 1631 B->ops->mult = MatMult_SeqSBAIJ_7; 1632 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1633 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1634 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1635 break; 1636 } 1637 } 1638 1639 b->mbs = mbs; 1640 b->nbs = nbs; 1641 if (!skipallocation) { 1642 if (!b->imax) { 1643 ierr = PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);CHKERRQ(ierr); 1644 1645 b->free_imax_ilen = PETSC_TRUE; 1646 1647 ierr = PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1648 } 1649 if (!nnz) { 1650 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1651 else if (nz <= 0) nz = 1; 1652 for (i=0; i<mbs; i++) b->imax[i] = nz; 1653 nz = nz*mbs; /* total nz */ 1654 } else { 1655 nz = 0; 1656 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1657 } 1658 /* b->ilen will count nonzeros in each block row so far. */ 1659 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1660 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1661 1662 /* allocate the matrix space */ 1663 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1664 ierr = PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);CHKERRQ(ierr); 1665 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1666 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1667 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1668 1669 b->singlemalloc = PETSC_TRUE; 1670 1671 /* pointer to beginning of each row */ 1672 b->i[0] = 0; 1673 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1674 1675 b->free_a = PETSC_TRUE; 1676 b->free_ij = PETSC_TRUE; 1677 } else { 1678 b->free_a = PETSC_FALSE; 1679 b->free_ij = PETSC_FALSE; 1680 } 1681 1682 B->rmap->bs = bs; 1683 b->bs2 = bs2; 1684 b->nz = 0; 1685 b->maxnz = nz; 1686 1687 b->inew = 0; 1688 b->jnew = 0; 1689 b->anew = 0; 1690 b->a2anew = 0; 1691 b->permute = PETSC_FALSE; 1692 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1693 PetscFunctionReturn(0); 1694 } 1695 1696 #undef __FUNCT__ 1697 #define __FUNCT__ "MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ" 1698 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1699 { 1700 PetscInt i,j,m,nz,nz_max=0,*nnz; 1701 PetscScalar *values=0; 1702 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1703 PetscErrorCode ierr; 1704 PetscFunctionBegin; 1705 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs); 1706 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1707 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1708 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1709 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1710 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1711 m = B->rmap->n/bs; 1712 1713 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]); 1714 ierr = PetscMalloc1(m+1,&nnz);CHKERRQ(ierr); 1715 for (i=0; i<m; i++) { 1716 nz = ii[i+1] - ii[i]; 1717 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz); 1718 nz_max = PetscMax(nz_max,nz); 1719 nnz[i] = nz; 1720 } 1721 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1722 ierr = PetscFree(nnz);CHKERRQ(ierr); 1723 1724 values = (PetscScalar*)V; 1725 if (!values) { 1726 ierr = PetscCalloc1(bs*bs*nz_max,&values);CHKERRQ(ierr); 1727 } 1728 for (i=0; i<m; i++) { 1729 PetscInt ncols = ii[i+1] - ii[i]; 1730 const PetscInt *icols = jj + ii[i]; 1731 if (!roworiented || bs == 1) { 1732 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1733 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr); 1734 } else { 1735 for (j=0; j<ncols; j++) { 1736 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1737 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr); 1738 } 1739 } 1740 } 1741 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); } 1742 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1743 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1744 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1745 PetscFunctionReturn(0); 1746 } 1747 1748 /* 1749 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1750 */ 1751 #undef __FUNCT__ 1752 #define __FUNCT__ "MatSeqSBAIJSetNumericFactorization_inplace" 1753 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1754 { 1755 PetscErrorCode ierr; 1756 PetscBool flg = PETSC_FALSE; 1757 PetscInt bs = B->rmap->bs; 1758 1759 PetscFunctionBegin; 1760 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1761 if (flg) bs = 8; 1762 1763 if (!natural) { 1764 switch (bs) { 1765 case 1: 1766 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1767 break; 1768 case 2: 1769 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1770 break; 1771 case 3: 1772 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1773 break; 1774 case 4: 1775 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1776 break; 1777 case 5: 1778 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1779 break; 1780 case 6: 1781 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1782 break; 1783 case 7: 1784 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1785 break; 1786 default: 1787 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1788 break; 1789 } 1790 } else { 1791 switch (bs) { 1792 case 1: 1793 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1794 break; 1795 case 2: 1796 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1797 break; 1798 case 3: 1799 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1800 break; 1801 case 4: 1802 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1803 break; 1804 case 5: 1805 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1806 break; 1807 case 6: 1808 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1809 break; 1810 case 7: 1811 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1812 break; 1813 default: 1814 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1815 break; 1816 } 1817 } 1818 PetscFunctionReturn(0); 1819 } 1820 1821 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1822 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1823 1824 #undef __FUNCT__ 1825 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1826 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1827 { 1828 PetscInt n = A->rmap->n; 1829 PetscErrorCode ierr; 1830 1831 PetscFunctionBegin; 1832 #if defined(PETSC_USE_COMPLEX) 1833 if (A->hermitian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported"); 1834 #endif 1835 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 1836 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1837 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1838 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1839 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1840 1841 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1842 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1843 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1844 1845 (*B)->factortype = ftype; 1846 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 1847 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 1848 PetscFunctionReturn(0); 1849 } 1850 1851 /*MC 1852 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1853 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1854 1855 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1856 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1857 1858 Options Database Keys: 1859 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1860 1861 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1862 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1863 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. 1864 1865 1866 Level: beginner 1867 1868 .seealso: MatCreateSeqSBAIJ 1869 M*/ 1870 1871 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1872 1873 #undef __FUNCT__ 1874 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1875 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1876 { 1877 Mat_SeqSBAIJ *b; 1878 PetscErrorCode ierr; 1879 PetscMPIInt size; 1880 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1881 1882 PetscFunctionBegin; 1883 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 1884 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1885 1886 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1887 B->data = (void*)b; 1888 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1889 1890 B->ops->destroy = MatDestroy_SeqSBAIJ; 1891 B->ops->view = MatView_SeqSBAIJ; 1892 b->row = 0; 1893 b->icol = 0; 1894 b->reallocs = 0; 1895 b->saved_values = 0; 1896 b->inode.limit = 5; 1897 b->inode.max_limit = 5; 1898 1899 b->roworiented = PETSC_TRUE; 1900 b->nonew = 0; 1901 b->diag = 0; 1902 b->solve_work = 0; 1903 b->mult_work = 0; 1904 B->spptr = 0; 1905 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1906 b->keepnonzeropattern = PETSC_FALSE; 1907 1908 b->inew = 0; 1909 b->jnew = 0; 1910 b->anew = 0; 1911 b->a2anew = 0; 1912 b->permute = PETSC_FALSE; 1913 1914 b->ignore_ltriangular = PETSC_TRUE; 1915 1916 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1917 1918 b->getrow_utriangular = PETSC_FALSE; 1919 1920 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 1921 1922 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1923 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1924 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1925 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1926 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1927 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1928 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr); 1929 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqsbstrm_C",MatConvert_SeqSBAIJ_SeqSBSTRM);CHKERRQ(ierr); 1930 #if defined(PETSC_HAVE_ELEMENTAL) 1931 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);CHKERRQ(ierr); 1932 #endif 1933 1934 B->symmetric = PETSC_TRUE; 1935 B->structurally_symmetric = PETSC_TRUE; 1936 B->symmetric_set = PETSC_TRUE; 1937 B->structurally_symmetric_set = PETSC_TRUE; 1938 1939 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1940 1941 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1942 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 1943 if (no_unroll) { 1944 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 1945 } 1946 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 1947 if (no_inode) { 1948 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 1949 } 1950 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 1951 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1952 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1953 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1954 PetscFunctionReturn(0); 1955 } 1956 1957 #undef __FUNCT__ 1958 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1959 /*@C 1960 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1961 compressed row) format. For good matrix assembly performance the 1962 user should preallocate the matrix storage by setting the parameter nz 1963 (or the array nnz). By setting these parameters accurately, performance 1964 during matrix assembly can be increased by more than a factor of 50. 1965 1966 Collective on Mat 1967 1968 Input Parameters: 1969 + B - the symmetric matrix 1970 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 1971 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 1972 . nz - number of block nonzeros per block row (same for all rows) 1973 - nnz - array containing the number of block nonzeros in the upper triangular plus 1974 diagonal portion of each block (possibly different for each block row) or NULL 1975 1976 Options Database Keys: 1977 . -mat_no_unroll - uses code that does not unroll the loops in the 1978 block calculations (much slower) 1979 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1980 1981 Level: intermediate 1982 1983 Notes: 1984 Specify the preallocated storage with either nz or nnz (not both). 1985 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 1986 allocation. See Users-Manual: ch_mat for details. 1987 1988 You can call MatGetInfo() to get information on how effective the preallocation was; 1989 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1990 You can also run with the option -info and look for messages with the string 1991 malloc in them to see if additional memory allocation was needed. 1992 1993 If the nnz parameter is given then the nz parameter is ignored 1994 1995 1996 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 1997 @*/ 1998 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1999 { 2000 PetscErrorCode ierr; 2001 2002 PetscFunctionBegin; 2003 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2004 PetscValidType(B,1); 2005 PetscValidLogicalCollectiveInt(B,bs,2); 2006 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2007 PetscFunctionReturn(0); 2008 } 2009 2010 #undef __FUNCT__ 2011 #define __FUNCT__ "MatSeqSBAIJSetPreallocationCSR" 2012 /*@C 2013 MatSeqSBAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in symmetric block AIJ format. 2014 2015 Input Parameters: 2016 + B - the matrix 2017 . bs - size of block, the blocks are ALWAYS square. 2018 . i - the indices into j for the start of each local row (starts with zero) 2019 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2020 - v - optional values in the matrix 2021 2022 Level: developer 2023 2024 Notes: 2025 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2026 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2027 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2028 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2029 block column and the second index is over columns within a block. 2030 2031 .keywords: matrix, block, aij, compressed row, sparse 2032 2033 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ 2034 @*/ 2035 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2036 { 2037 PetscErrorCode ierr; 2038 2039 PetscFunctionBegin; 2040 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2041 PetscValidType(B,1); 2042 PetscValidLogicalCollectiveInt(B,bs,2); 2043 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2044 PetscFunctionReturn(0); 2045 } 2046 2047 #undef __FUNCT__ 2048 #define __FUNCT__ "MatCreateSeqSBAIJ" 2049 /*@C 2050 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2051 compressed row) format. For good matrix assembly performance the 2052 user should preallocate the matrix storage by setting the parameter nz 2053 (or the array nnz). By setting these parameters accurately, performance 2054 during matrix assembly can be increased by more than a factor of 50. 2055 2056 Collective on MPI_Comm 2057 2058 Input Parameters: 2059 + comm - MPI communicator, set to PETSC_COMM_SELF 2060 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2061 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2062 . m - number of rows, or number of columns 2063 . nz - number of block nonzeros per block row (same for all rows) 2064 - nnz - array containing the number of block nonzeros in the upper triangular plus 2065 diagonal portion of each block (possibly different for each block row) or NULL 2066 2067 Output Parameter: 2068 . A - the symmetric matrix 2069 2070 Options Database Keys: 2071 . -mat_no_unroll - uses code that does not unroll the loops in the 2072 block calculations (much slower) 2073 . -mat_block_size - size of the blocks to use 2074 2075 Level: intermediate 2076 2077 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2078 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2079 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2080 2081 Notes: 2082 The number of rows and columns must be divisible by blocksize. 2083 This matrix type does not support complex Hermitian operation. 2084 2085 Specify the preallocated storage with either nz or nnz (not both). 2086 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2087 allocation. See Users-Manual: ch_mat for details. 2088 2089 If the nnz parameter is given then the nz parameter is ignored 2090 2091 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2092 @*/ 2093 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2094 { 2095 PetscErrorCode ierr; 2096 2097 PetscFunctionBegin; 2098 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2099 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2100 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2101 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2102 PetscFunctionReturn(0); 2103 } 2104 2105 #undef __FUNCT__ 2106 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2107 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2108 { 2109 Mat C; 2110 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2111 PetscErrorCode ierr; 2112 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2113 2114 PetscFunctionBegin; 2115 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2116 2117 *B = 0; 2118 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2119 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2120 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2121 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2122 c = (Mat_SeqSBAIJ*)C->data; 2123 2124 C->preallocated = PETSC_TRUE; 2125 C->factortype = A->factortype; 2126 c->row = 0; 2127 c->icol = 0; 2128 c->saved_values = 0; 2129 c->keepnonzeropattern = a->keepnonzeropattern; 2130 C->assembled = PETSC_TRUE; 2131 2132 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2133 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2134 c->bs2 = a->bs2; 2135 c->mbs = a->mbs; 2136 c->nbs = a->nbs; 2137 2138 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2139 c->imax = a->imax; 2140 c->ilen = a->ilen; 2141 c->free_imax_ilen = PETSC_FALSE; 2142 } else { 2143 ierr = PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);CHKERRQ(ierr); 2144 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2145 for (i=0; i<mbs; i++) { 2146 c->imax[i] = a->imax[i]; 2147 c->ilen[i] = a->ilen[i]; 2148 } 2149 c->free_imax_ilen = PETSC_TRUE; 2150 } 2151 2152 /* allocate the matrix space */ 2153 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2154 ierr = PetscMalloc1(bs2*nz,&c->a);CHKERRQ(ierr); 2155 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2156 c->i = a->i; 2157 c->j = a->j; 2158 c->singlemalloc = PETSC_FALSE; 2159 c->free_a = PETSC_TRUE; 2160 c->free_ij = PETSC_FALSE; 2161 c->parent = A; 2162 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2163 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2164 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2165 } else { 2166 ierr = PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);CHKERRQ(ierr); 2167 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2168 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2169 c->singlemalloc = PETSC_TRUE; 2170 c->free_a = PETSC_TRUE; 2171 c->free_ij = PETSC_TRUE; 2172 } 2173 if (mbs > 0) { 2174 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2175 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2176 } 2177 if (cpvalues == MAT_COPY_VALUES) { 2178 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2179 } else { 2180 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2181 } 2182 if (a->jshort) { 2183 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2184 /* if the parent matrix is reassembled, this child matrix will never notice */ 2185 ierr = PetscMalloc1(nz,&c->jshort);CHKERRQ(ierr); 2186 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2187 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2188 2189 c->free_jshort = PETSC_TRUE; 2190 } 2191 } 2192 2193 c->roworiented = a->roworiented; 2194 c->nonew = a->nonew; 2195 2196 if (a->diag) { 2197 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2198 c->diag = a->diag; 2199 c->free_diag = PETSC_FALSE; 2200 } else { 2201 ierr = PetscMalloc1(mbs,&c->diag);CHKERRQ(ierr); 2202 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2203 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2204 c->free_diag = PETSC_TRUE; 2205 } 2206 } 2207 c->nz = a->nz; 2208 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2209 c->solve_work = 0; 2210 c->mult_work = 0; 2211 2212 *B = C; 2213 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2214 PetscFunctionReturn(0); 2215 } 2216 2217 #undef __FUNCT__ 2218 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2219 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2220 { 2221 Mat_SeqSBAIJ *a; 2222 PetscErrorCode ierr; 2223 int fd; 2224 PetscMPIInt size; 2225 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs = newmat->rmap->bs; 2226 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2227 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2228 PetscInt *masked,nmask,tmp,bs2,ishift; 2229 PetscScalar *aa; 2230 MPI_Comm comm; 2231 2232 PetscFunctionBegin; 2233 /* force binary viewer to load .info file if it has not yet done so */ 2234 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 2235 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2236 ierr = PetscOptionsGetInt(((PetscObject)newmat)->options,((PetscObject)newmat)->prefix,"-matload_block_size",&bs,NULL);CHKERRQ(ierr); 2237 if (bs < 0) bs = 1; 2238 bs2 = bs*bs; 2239 2240 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2241 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2242 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2243 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2244 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2245 M = header[1]; N = header[2]; nz = header[3]; 2246 2247 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2248 2249 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2250 2251 /* 2252 This code adds extra rows to make sure the number of rows is 2253 divisible by the blocksize 2254 */ 2255 mbs = M/bs; 2256 extra_rows = bs - M + bs*(mbs); 2257 if (extra_rows == bs) extra_rows = 0; 2258 else mbs++; 2259 if (extra_rows) { 2260 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2261 } 2262 2263 /* Set global sizes if not already set */ 2264 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2265 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2266 } else { /* Check if the matrix global sizes are correct */ 2267 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2268 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); 2269 } 2270 2271 /* read in row lengths */ 2272 ierr = PetscMalloc1(M+extra_rows,&rowlengths);CHKERRQ(ierr); 2273 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2274 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2275 2276 /* read in column indices */ 2277 ierr = PetscMalloc1(nz+extra_rows,&jj);CHKERRQ(ierr); 2278 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2279 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2280 2281 /* loop over row lengths determining block row lengths */ 2282 ierr = PetscCalloc1(mbs,&s_browlengths);CHKERRQ(ierr); 2283 ierr = PetscMalloc2(mbs,&mask,mbs,&masked);CHKERRQ(ierr); 2284 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2285 rowcount = 0; 2286 nzcount = 0; 2287 for (i=0; i<mbs; i++) { 2288 nmask = 0; 2289 for (j=0; j<bs; j++) { 2290 kmax = rowlengths[rowcount]; 2291 for (k=0; k<kmax; k++) { 2292 tmp = jj[nzcount++]/bs; /* block col. index */ 2293 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2294 } 2295 rowcount++; 2296 } 2297 s_browlengths[i] += nmask; 2298 2299 /* zero out the mask elements we set */ 2300 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2301 } 2302 2303 /* Do preallocation */ 2304 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2305 a = (Mat_SeqSBAIJ*)newmat->data; 2306 2307 /* set matrix "i" values */ 2308 a->i[0] = 0; 2309 for (i=1; i<= mbs; i++) { 2310 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2311 a->ilen[i-1] = s_browlengths[i-1]; 2312 } 2313 a->nz = a->i[mbs]; 2314 2315 /* read in nonzero values */ 2316 ierr = PetscMalloc1(nz+extra_rows,&aa);CHKERRQ(ierr); 2317 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2318 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2319 2320 /* set "a" and "j" values into matrix */ 2321 nzcount = 0; jcount = 0; 2322 for (i=0; i<mbs; i++) { 2323 nzcountb = nzcount; 2324 nmask = 0; 2325 for (j=0; j<bs; j++) { 2326 kmax = rowlengths[i*bs+j]; 2327 for (k=0; k<kmax; k++) { 2328 tmp = jj[nzcount++]/bs; /* block col. index */ 2329 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2330 } 2331 } 2332 /* sort the masked values */ 2333 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2334 2335 /* set "j" values into matrix */ 2336 maskcount = 1; 2337 for (j=0; j<nmask; j++) { 2338 a->j[jcount++] = masked[j]; 2339 mask[masked[j]] = maskcount++; 2340 } 2341 2342 /* set "a" values into matrix */ 2343 ishift = bs2*a->i[i]; 2344 for (j=0; j<bs; j++) { 2345 kmax = rowlengths[i*bs+j]; 2346 for (k=0; k<kmax; k++) { 2347 tmp = jj[nzcountb]/bs; /* block col. index */ 2348 if (tmp >= i) { 2349 block = mask[tmp] - 1; 2350 point = jj[nzcountb] - bs*tmp; 2351 idx = ishift + bs2*block + j + bs*point; 2352 a->a[idx] = aa[nzcountb]; 2353 } 2354 nzcountb++; 2355 } 2356 } 2357 /* zero out the mask elements we set */ 2358 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2359 } 2360 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2361 2362 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2363 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2364 ierr = PetscFree(aa);CHKERRQ(ierr); 2365 ierr = PetscFree(jj);CHKERRQ(ierr); 2366 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2367 2368 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2369 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2370 PetscFunctionReturn(0); 2371 } 2372 2373 #undef __FUNCT__ 2374 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2375 /*@ 2376 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2377 (upper triangular entries in CSR format) provided by the user. 2378 2379 Collective on MPI_Comm 2380 2381 Input Parameters: 2382 + comm - must be an MPI communicator of size 1 2383 . bs - size of block 2384 . m - number of rows 2385 . n - number of columns 2386 . i - row indices 2387 . j - column indices 2388 - a - matrix values 2389 2390 Output Parameter: 2391 . mat - the matrix 2392 2393 Level: advanced 2394 2395 Notes: 2396 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2397 once the matrix is destroyed 2398 2399 You cannot set new nonzero locations into this matrix, that will generate an error. 2400 2401 The i and j indices are 0 based 2402 2403 When block size is greater than 1 the matrix values must be stored using the SBAIJ storage format (see the SBAIJ code to determine this). For block size of 1 2404 it is the regular CSR format excluding the lower triangular elements. 2405 2406 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2407 2408 @*/ 2409 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 2410 { 2411 PetscErrorCode ierr; 2412 PetscInt ii; 2413 Mat_SeqSBAIJ *sbaij; 2414 2415 PetscFunctionBegin; 2416 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2417 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2418 2419 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2420 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2421 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2422 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2423 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2424 ierr = PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);CHKERRQ(ierr); 2425 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2426 2427 sbaij->i = i; 2428 sbaij->j = j; 2429 sbaij->a = a; 2430 2431 sbaij->singlemalloc = PETSC_FALSE; 2432 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2433 sbaij->free_a = PETSC_FALSE; 2434 sbaij->free_ij = PETSC_FALSE; 2435 2436 for (ii=0; ii<m; ii++) { 2437 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2438 #if defined(PETSC_USE_DEBUG) 2439 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2440 #endif 2441 } 2442 #if defined(PETSC_USE_DEBUG) 2443 for (ii=0; ii<sbaij->i[m]; ii++) { 2444 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2445 if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2446 } 2447 #endif 2448 2449 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2450 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2451 PetscFunctionReturn(0); 2452 } 2453 2454 #undef __FUNCT__ 2455 #define __FUNCT__ "MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ" 2456 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2457 { 2458 PetscErrorCode ierr; 2459 2460 PetscFunctionBegin; 2461 ierr = MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2462 PetscFunctionReturn(0); 2463 } 2464 2465 2466 2467 2468