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 (*B)->factortype = ftype; 1845 PetscFunctionReturn(0); 1846 } 1847 1848 /*MC 1849 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1850 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1851 1852 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1853 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1854 1855 Options Database Keys: 1856 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1857 1858 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1859 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1860 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. 1861 1862 1863 Level: beginner 1864 1865 .seealso: MatCreateSeqSBAIJ 1866 M*/ 1867 1868 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1869 1870 #undef __FUNCT__ 1871 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1872 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1873 { 1874 Mat_SeqSBAIJ *b; 1875 PetscErrorCode ierr; 1876 PetscMPIInt size; 1877 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1878 1879 PetscFunctionBegin; 1880 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 1881 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1882 1883 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1884 B->data = (void*)b; 1885 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1886 1887 B->ops->destroy = MatDestroy_SeqSBAIJ; 1888 B->ops->view = MatView_SeqSBAIJ; 1889 b->row = 0; 1890 b->icol = 0; 1891 b->reallocs = 0; 1892 b->saved_values = 0; 1893 b->inode.limit = 5; 1894 b->inode.max_limit = 5; 1895 1896 b->roworiented = PETSC_TRUE; 1897 b->nonew = 0; 1898 b->diag = 0; 1899 b->solve_work = 0; 1900 b->mult_work = 0; 1901 B->spptr = 0; 1902 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1903 b->keepnonzeropattern = PETSC_FALSE; 1904 1905 b->inew = 0; 1906 b->jnew = 0; 1907 b->anew = 0; 1908 b->a2anew = 0; 1909 b->permute = PETSC_FALSE; 1910 1911 b->ignore_ltriangular = PETSC_TRUE; 1912 1913 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1914 1915 b->getrow_utriangular = PETSC_FALSE; 1916 1917 ierr = PetscOptionsGetBool(((PetscObject)B)->options,((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 1918 1919 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1920 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1921 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1922 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1923 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1924 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1925 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr); 1926 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqsbstrm_C",MatConvert_SeqSBAIJ_SeqSBSTRM);CHKERRQ(ierr); 1927 #if defined(PETSC_HAVE_ELEMENTAL) 1928 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);CHKERRQ(ierr); 1929 #endif 1930 1931 B->symmetric = PETSC_TRUE; 1932 B->structurally_symmetric = PETSC_TRUE; 1933 B->symmetric_set = PETSC_TRUE; 1934 B->structurally_symmetric_set = PETSC_TRUE; 1935 1936 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1937 1938 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1939 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 1940 if (no_unroll) { 1941 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 1942 } 1943 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 1944 if (no_inode) { 1945 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 1946 } 1947 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 1948 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1949 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1950 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1951 PetscFunctionReturn(0); 1952 } 1953 1954 #undef __FUNCT__ 1955 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1956 /*@C 1957 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1958 compressed row) format. For good matrix assembly performance the 1959 user should preallocate the matrix storage by setting the parameter nz 1960 (or the array nnz). By setting these parameters accurately, performance 1961 during matrix assembly can be increased by more than a factor of 50. 1962 1963 Collective on Mat 1964 1965 Input Parameters: 1966 + B - the symmetric matrix 1967 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 1968 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 1969 . nz - number of block nonzeros per block row (same for all rows) 1970 - nnz - array containing the number of block nonzeros in the upper triangular plus 1971 diagonal portion of each block (possibly different for each block row) or NULL 1972 1973 Options Database Keys: 1974 . -mat_no_unroll - uses code that does not unroll the loops in the 1975 block calculations (much slower) 1976 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1977 1978 Level: intermediate 1979 1980 Notes: 1981 Specify the preallocated storage with either nz or nnz (not both). 1982 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 1983 allocation. See Users-Manual: ch_mat for details. 1984 1985 You can call MatGetInfo() to get information on how effective the preallocation was; 1986 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1987 You can also run with the option -info and look for messages with the string 1988 malloc in them to see if additional memory allocation was needed. 1989 1990 If the nnz parameter is given then the nz parameter is ignored 1991 1992 1993 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 1994 @*/ 1995 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1996 { 1997 PetscErrorCode ierr; 1998 1999 PetscFunctionBegin; 2000 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2001 PetscValidType(B,1); 2002 PetscValidLogicalCollectiveInt(B,bs,2); 2003 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2004 PetscFunctionReturn(0); 2005 } 2006 2007 #undef __FUNCT__ 2008 #define __FUNCT__ "MatSeqSBAIJSetPreallocationCSR" 2009 /*@C 2010 MatSeqSBAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in symmetric block AIJ format. 2011 2012 Input Parameters: 2013 + B - the matrix 2014 . bs - size of block, the blocks are ALWAYS square. 2015 . i - the indices into j for the start of each local row (starts with zero) 2016 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2017 - v - optional values in the matrix 2018 2019 Level: developer 2020 2021 Notes: 2022 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2023 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2024 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2025 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2026 block column and the second index is over columns within a block. 2027 2028 .keywords: matrix, block, aij, compressed row, sparse 2029 2030 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ 2031 @*/ 2032 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2033 { 2034 PetscErrorCode ierr; 2035 2036 PetscFunctionBegin; 2037 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2038 PetscValidType(B,1); 2039 PetscValidLogicalCollectiveInt(B,bs,2); 2040 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2041 PetscFunctionReturn(0); 2042 } 2043 2044 #undef __FUNCT__ 2045 #define __FUNCT__ "MatCreateSeqSBAIJ" 2046 /*@C 2047 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2048 compressed row) format. For good matrix assembly performance the 2049 user should preallocate the matrix storage by setting the parameter nz 2050 (or the array nnz). By setting these parameters accurately, performance 2051 during matrix assembly can be increased by more than a factor of 50. 2052 2053 Collective on MPI_Comm 2054 2055 Input Parameters: 2056 + comm - MPI communicator, set to PETSC_COMM_SELF 2057 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2058 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2059 . m - number of rows, or number of columns 2060 . nz - number of block nonzeros per block row (same for all rows) 2061 - nnz - array containing the number of block nonzeros in the upper triangular plus 2062 diagonal portion of each block (possibly different for each block row) or NULL 2063 2064 Output Parameter: 2065 . A - the symmetric matrix 2066 2067 Options Database Keys: 2068 . -mat_no_unroll - uses code that does not unroll the loops in the 2069 block calculations (much slower) 2070 . -mat_block_size - size of the blocks to use 2071 2072 Level: intermediate 2073 2074 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2075 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2076 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2077 2078 Notes: 2079 The number of rows and columns must be divisible by blocksize. 2080 This matrix type does not support complex Hermitian operation. 2081 2082 Specify the preallocated storage with either nz or nnz (not both). 2083 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2084 allocation. See Users-Manual: ch_mat for details. 2085 2086 If the nnz parameter is given then the nz parameter is ignored 2087 2088 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2089 @*/ 2090 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2091 { 2092 PetscErrorCode ierr; 2093 2094 PetscFunctionBegin; 2095 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2096 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2097 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2098 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2099 PetscFunctionReturn(0); 2100 } 2101 2102 #undef __FUNCT__ 2103 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2104 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2105 { 2106 Mat C; 2107 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2108 PetscErrorCode ierr; 2109 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2110 2111 PetscFunctionBegin; 2112 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2113 2114 *B = 0; 2115 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2116 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2117 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2118 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2119 c = (Mat_SeqSBAIJ*)C->data; 2120 2121 C->preallocated = PETSC_TRUE; 2122 C->factortype = A->factortype; 2123 c->row = 0; 2124 c->icol = 0; 2125 c->saved_values = 0; 2126 c->keepnonzeropattern = a->keepnonzeropattern; 2127 C->assembled = PETSC_TRUE; 2128 2129 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2130 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2131 c->bs2 = a->bs2; 2132 c->mbs = a->mbs; 2133 c->nbs = a->nbs; 2134 2135 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2136 c->imax = a->imax; 2137 c->ilen = a->ilen; 2138 c->free_imax_ilen = PETSC_FALSE; 2139 } else { 2140 ierr = PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);CHKERRQ(ierr); 2141 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2142 for (i=0; i<mbs; i++) { 2143 c->imax[i] = a->imax[i]; 2144 c->ilen[i] = a->ilen[i]; 2145 } 2146 c->free_imax_ilen = PETSC_TRUE; 2147 } 2148 2149 /* allocate the matrix space */ 2150 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2151 ierr = PetscMalloc1(bs2*nz,&c->a);CHKERRQ(ierr); 2152 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2153 c->i = a->i; 2154 c->j = a->j; 2155 c->singlemalloc = PETSC_FALSE; 2156 c->free_a = PETSC_TRUE; 2157 c->free_ij = PETSC_FALSE; 2158 c->parent = A; 2159 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2160 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2161 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2162 } else { 2163 ierr = PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);CHKERRQ(ierr); 2164 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2165 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2166 c->singlemalloc = PETSC_TRUE; 2167 c->free_a = PETSC_TRUE; 2168 c->free_ij = PETSC_TRUE; 2169 } 2170 if (mbs > 0) { 2171 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2172 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2173 } 2174 if (cpvalues == MAT_COPY_VALUES) { 2175 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2176 } else { 2177 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2178 } 2179 if (a->jshort) { 2180 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2181 /* if the parent matrix is reassembled, this child matrix will never notice */ 2182 ierr = PetscMalloc1(nz,&c->jshort);CHKERRQ(ierr); 2183 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2184 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2185 2186 c->free_jshort = PETSC_TRUE; 2187 } 2188 } 2189 2190 c->roworiented = a->roworiented; 2191 c->nonew = a->nonew; 2192 2193 if (a->diag) { 2194 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2195 c->diag = a->diag; 2196 c->free_diag = PETSC_FALSE; 2197 } else { 2198 ierr = PetscMalloc1(mbs,&c->diag);CHKERRQ(ierr); 2199 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2200 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2201 c->free_diag = PETSC_TRUE; 2202 } 2203 } 2204 c->nz = a->nz; 2205 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2206 c->solve_work = 0; 2207 c->mult_work = 0; 2208 2209 *B = C; 2210 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2211 PetscFunctionReturn(0); 2212 } 2213 2214 #undef __FUNCT__ 2215 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2216 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2217 { 2218 Mat_SeqSBAIJ *a; 2219 PetscErrorCode ierr; 2220 int fd; 2221 PetscMPIInt size; 2222 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs = newmat->rmap->bs; 2223 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2224 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2225 PetscInt *masked,nmask,tmp,bs2,ishift; 2226 PetscScalar *aa; 2227 MPI_Comm comm; 2228 2229 PetscFunctionBegin; 2230 /* force binary viewer to load .info file if it has not yet done so */ 2231 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 2232 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2233 ierr = PetscOptionsGetInt(((PetscObject)newmat)->options,((PetscObject)newmat)->prefix,"-matload_block_size",&bs,NULL);CHKERRQ(ierr); 2234 if (bs < 0) bs = 1; 2235 bs2 = bs*bs; 2236 2237 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2238 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2239 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2240 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2241 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2242 M = header[1]; N = header[2]; nz = header[3]; 2243 2244 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2245 2246 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2247 2248 /* 2249 This code adds extra rows to make sure the number of rows is 2250 divisible by the blocksize 2251 */ 2252 mbs = M/bs; 2253 extra_rows = bs - M + bs*(mbs); 2254 if (extra_rows == bs) extra_rows = 0; 2255 else mbs++; 2256 if (extra_rows) { 2257 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2258 } 2259 2260 /* Set global sizes if not already set */ 2261 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2262 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2263 } else { /* Check if the matrix global sizes are correct */ 2264 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2265 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); 2266 } 2267 2268 /* read in row lengths */ 2269 ierr = PetscMalloc1(M+extra_rows,&rowlengths);CHKERRQ(ierr); 2270 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2271 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2272 2273 /* read in column indices */ 2274 ierr = PetscMalloc1(nz+extra_rows,&jj);CHKERRQ(ierr); 2275 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2276 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2277 2278 /* loop over row lengths determining block row lengths */ 2279 ierr = PetscCalloc1(mbs,&s_browlengths);CHKERRQ(ierr); 2280 ierr = PetscMalloc2(mbs,&mask,mbs,&masked);CHKERRQ(ierr); 2281 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2282 rowcount = 0; 2283 nzcount = 0; 2284 for (i=0; i<mbs; i++) { 2285 nmask = 0; 2286 for (j=0; j<bs; j++) { 2287 kmax = rowlengths[rowcount]; 2288 for (k=0; k<kmax; k++) { 2289 tmp = jj[nzcount++]/bs; /* block col. index */ 2290 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2291 } 2292 rowcount++; 2293 } 2294 s_browlengths[i] += nmask; 2295 2296 /* zero out the mask elements we set */ 2297 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2298 } 2299 2300 /* Do preallocation */ 2301 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2302 a = (Mat_SeqSBAIJ*)newmat->data; 2303 2304 /* set matrix "i" values */ 2305 a->i[0] = 0; 2306 for (i=1; i<= mbs; i++) { 2307 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2308 a->ilen[i-1] = s_browlengths[i-1]; 2309 } 2310 a->nz = a->i[mbs]; 2311 2312 /* read in nonzero values */ 2313 ierr = PetscMalloc1(nz+extra_rows,&aa);CHKERRQ(ierr); 2314 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2315 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2316 2317 /* set "a" and "j" values into matrix */ 2318 nzcount = 0; jcount = 0; 2319 for (i=0; i<mbs; i++) { 2320 nzcountb = nzcount; 2321 nmask = 0; 2322 for (j=0; j<bs; j++) { 2323 kmax = rowlengths[i*bs+j]; 2324 for (k=0; k<kmax; k++) { 2325 tmp = jj[nzcount++]/bs; /* block col. index */ 2326 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2327 } 2328 } 2329 /* sort the masked values */ 2330 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2331 2332 /* set "j" values into matrix */ 2333 maskcount = 1; 2334 for (j=0; j<nmask; j++) { 2335 a->j[jcount++] = masked[j]; 2336 mask[masked[j]] = maskcount++; 2337 } 2338 2339 /* set "a" values into matrix */ 2340 ishift = bs2*a->i[i]; 2341 for (j=0; j<bs; j++) { 2342 kmax = rowlengths[i*bs+j]; 2343 for (k=0; k<kmax; k++) { 2344 tmp = jj[nzcountb]/bs; /* block col. index */ 2345 if (tmp >= i) { 2346 block = mask[tmp] - 1; 2347 point = jj[nzcountb] - bs*tmp; 2348 idx = ishift + bs2*block + j + bs*point; 2349 a->a[idx] = aa[nzcountb]; 2350 } 2351 nzcountb++; 2352 } 2353 } 2354 /* zero out the mask elements we set */ 2355 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2356 } 2357 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2358 2359 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2360 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2361 ierr = PetscFree(aa);CHKERRQ(ierr); 2362 ierr = PetscFree(jj);CHKERRQ(ierr); 2363 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2364 2365 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2366 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2367 PetscFunctionReturn(0); 2368 } 2369 2370 #undef __FUNCT__ 2371 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2372 /*@ 2373 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2374 (upper triangular entries in CSR format) provided by the user. 2375 2376 Collective on MPI_Comm 2377 2378 Input Parameters: 2379 + comm - must be an MPI communicator of size 1 2380 . bs - size of block 2381 . m - number of rows 2382 . n - number of columns 2383 . i - row indices 2384 . j - column indices 2385 - a - matrix values 2386 2387 Output Parameter: 2388 . mat - the matrix 2389 2390 Level: advanced 2391 2392 Notes: 2393 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2394 once the matrix is destroyed 2395 2396 You cannot set new nonzero locations into this matrix, that will generate an error. 2397 2398 The i and j indices are 0 based 2399 2400 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 2401 it is the regular CSR format excluding the lower triangular elements. 2402 2403 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2404 2405 @*/ 2406 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 2407 { 2408 PetscErrorCode ierr; 2409 PetscInt ii; 2410 Mat_SeqSBAIJ *sbaij; 2411 2412 PetscFunctionBegin; 2413 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2414 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2415 2416 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2417 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2418 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2419 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2420 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2421 ierr = PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);CHKERRQ(ierr); 2422 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2423 2424 sbaij->i = i; 2425 sbaij->j = j; 2426 sbaij->a = a; 2427 2428 sbaij->singlemalloc = PETSC_FALSE; 2429 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2430 sbaij->free_a = PETSC_FALSE; 2431 sbaij->free_ij = PETSC_FALSE; 2432 2433 for (ii=0; ii<m; ii++) { 2434 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2435 #if defined(PETSC_USE_DEBUG) 2436 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]); 2437 #endif 2438 } 2439 #if defined(PETSC_USE_DEBUG) 2440 for (ii=0; ii<sbaij->i[m]; ii++) { 2441 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2442 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]); 2443 } 2444 #endif 2445 2446 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2447 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2448 PetscFunctionReturn(0); 2449 } 2450 2451 #undef __FUNCT__ 2452 #define __FUNCT__ "MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ" 2453 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2454 { 2455 PetscErrorCode ierr; 2456 2457 PetscFunctionBegin; 2458 ierr = MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2459 PetscFunctionReturn(0); 2460 } 2461 2462 2463 2464 2465