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