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 if (A->factortype && bs>1) { 317 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); 318 PetscFunctionReturn(0); 319 } 320 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);CHKERRQ(ierr); 321 ierr = MatView(aij,viewer);CHKERRQ(ierr); 322 ierr = MatDestroy(&aij);CHKERRQ(ierr); 323 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 324 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 325 for (i=0; i<a->mbs; i++) { 326 for (j=0; j<bs; j++) { 327 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 328 for (k=a->i[i]; k<a->i[i+1]; k++) { 329 for (l=0; l<bs; l++) { 330 #if defined(PETSC_USE_COMPLEX) 331 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 332 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l, 333 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 334 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 335 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l, 336 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 337 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 338 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 339 } 340 #else 341 if (a->a[bs2*k + l*bs + j] != 0.0) { 342 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 343 } 344 #endif 345 } 346 } 347 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 348 } 349 } 350 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 351 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 352 PetscFunctionReturn(0); 353 } else { 354 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 355 if (A->factortype) { /* for factored matrix */ 356 if (bs>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet"); 357 358 diag=a->diag; 359 for (i=0; i<a->mbs; i++) { /* for row block i */ 360 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 361 /* diagonal entry */ 362 #if defined(PETSC_USE_COMPLEX) 363 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) { 364 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); 365 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) { 366 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); 367 } else { 368 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)PetscRealPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 369 } 370 #else 371 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[diag[i]],(double)(1.0/a->a[diag[i]]));CHKERRQ(ierr); 372 #endif 373 /* off-diagonal entries */ 374 for (k=a->i[i]; k<a->i[i+1]-1; k++) { 375 #if defined(PETSC_USE_COMPLEX) 376 if (PetscImaginaryPart(a->a[k]) > 0.0) { 377 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),(double)PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 378 } else if (PetscImaginaryPart(a->a[k]) < 0.0) { 379 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k],(double)PetscRealPart(a->a[k]),-(double)PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 380 } else { 381 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k],(double)PetscRealPart(a->a[k]));CHKERRQ(ierr); 382 } 383 #else 384 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[k],(double)a->a[k]);CHKERRQ(ierr); 385 #endif 386 } 387 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 388 } 389 390 } else { /* for non-factored matrix */ 391 for (i=0; i<a->mbs; i++) { /* for row block i */ 392 for (j=0; j<bs; j++) { /* for row bs*i + j */ 393 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 394 for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */ 395 for (l=0; l<bs; l++) { /* for column */ 396 #if defined(PETSC_USE_COMPLEX) 397 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 398 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",bs*a->j[k]+l, 399 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 400 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 401 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i) ",bs*a->j[k]+l, 402 (double)PetscRealPart(a->a[bs2*k + l*bs + j]),-(double)PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 403 } else { 404 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 405 } 406 #else 407 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",bs*a->j[k]+l,(double)a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 408 #endif 409 } 410 } 411 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 412 } 413 } 414 } 415 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 416 } 417 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 418 PetscFunctionReturn(0); 419 } 420 421 #include <petscdraw.h> 422 #undef __FUNCT__ 423 #define __FUNCT__ "MatView_SeqSBAIJ_Draw_Zoom" 424 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 425 { 426 Mat A = (Mat) Aa; 427 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 428 PetscErrorCode ierr; 429 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 430 PetscMPIInt rank; 431 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 432 MatScalar *aa; 433 MPI_Comm comm; 434 PetscViewer viewer; 435 436 PetscFunctionBegin; 437 /* 438 This is nasty. If this is called from an originally parallel matrix 439 then all processes call this,but only the first has the matrix so the 440 rest should return immediately. 441 */ 442 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 443 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 444 if (rank) PetscFunctionReturn(0); 445 446 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 447 448 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 449 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 450 451 /* loop over matrix elements drawing boxes */ 452 color = PETSC_DRAW_BLUE; 453 for (i=0,row=0; i<mbs; i++,row+=bs) { 454 for (j=a->i[i]; j<a->i[i+1]; j++) { 455 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 456 x_l = a->j[j]*bs; x_r = x_l + 1.0; 457 aa = a->a + j*bs2; 458 for (k=0; k<bs; k++) { 459 for (l=0; l<bs; l++) { 460 if (PetscRealPart(*aa++) >= 0.) continue; 461 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 462 } 463 } 464 } 465 } 466 color = PETSC_DRAW_CYAN; 467 for (i=0,row=0; i<mbs; i++,row+=bs) { 468 for (j=a->i[i]; j<a->i[i+1]; j++) { 469 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 470 x_l = a->j[j]*bs; x_r = x_l + 1.0; 471 aa = a->a + j*bs2; 472 for (k=0; k<bs; k++) { 473 for (l=0; l<bs; l++) { 474 if (PetscRealPart(*aa++) != 0.) continue; 475 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 476 } 477 } 478 } 479 } 480 481 color = PETSC_DRAW_RED; 482 for (i=0,row=0; i<mbs; i++,row+=bs) { 483 for (j=a->i[i]; j<a->i[i+1]; j++) { 484 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 485 x_l = a->j[j]*bs; x_r = x_l + 1.0; 486 aa = a->a + j*bs2; 487 for (k=0; k<bs; k++) { 488 for (l=0; l<bs; l++) { 489 if (PetscRealPart(*aa++) <= 0.) continue; 490 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 491 } 492 } 493 } 494 } 495 PetscFunctionReturn(0); 496 } 497 498 #undef __FUNCT__ 499 #define __FUNCT__ "MatView_SeqSBAIJ_Draw" 500 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 501 { 502 PetscErrorCode ierr; 503 PetscReal xl,yl,xr,yr,w,h; 504 PetscDraw draw; 505 PetscBool isnull; 506 507 PetscFunctionBegin; 508 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 509 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 510 511 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 512 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 513 xr += w; yr += h; xl = -w; yl = -h; 514 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 515 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 516 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 517 PetscFunctionReturn(0); 518 } 519 520 #undef __FUNCT__ 521 #define __FUNCT__ "MatView_SeqSBAIJ" 522 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 523 { 524 PetscErrorCode ierr; 525 PetscBool iascii,isdraw; 526 FILE *file = 0; 527 528 PetscFunctionBegin; 529 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 530 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 531 if (iascii) { 532 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 533 } else if (isdraw) { 534 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 535 } else { 536 Mat B; 537 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);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,"Row too large: row %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,"Column too large: col %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 nonzero (%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 /* -------------------------------------------------------------------*/ 1348 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1349 MatGetRow_SeqSBAIJ, 1350 MatRestoreRow_SeqSBAIJ, 1351 MatMult_SeqSBAIJ_N, 1352 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1353 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1354 MatMultAdd_SeqSBAIJ_N, 1355 0, 1356 0, 1357 0, 1358 /* 10*/ 0, 1359 0, 1360 MatCholeskyFactor_SeqSBAIJ, 1361 MatSOR_SeqSBAIJ, 1362 MatTranspose_SeqSBAIJ, 1363 /* 15*/ MatGetInfo_SeqSBAIJ, 1364 MatEqual_SeqSBAIJ, 1365 MatGetDiagonal_SeqSBAIJ, 1366 MatDiagonalScale_SeqSBAIJ, 1367 MatNorm_SeqSBAIJ, 1368 /* 20*/ 0, 1369 MatAssemblyEnd_SeqSBAIJ, 1370 MatSetOption_SeqSBAIJ, 1371 MatZeroEntries_SeqSBAIJ, 1372 /* 24*/ 0, 1373 0, 1374 0, 1375 0, 1376 0, 1377 /* 29*/ MatSetUp_SeqSBAIJ, 1378 0, 1379 0, 1380 0, 1381 0, 1382 /* 34*/ MatDuplicate_SeqSBAIJ, 1383 0, 1384 0, 1385 0, 1386 MatICCFactor_SeqSBAIJ, 1387 /* 39*/ MatAXPY_SeqSBAIJ, 1388 MatGetSubMatrices_SeqSBAIJ, 1389 MatIncreaseOverlap_SeqSBAIJ, 1390 MatGetValues_SeqSBAIJ, 1391 MatCopy_SeqSBAIJ, 1392 /* 44*/ 0, 1393 MatScale_SeqSBAIJ, 1394 0, 1395 0, 1396 MatZeroRowsColumns_SeqSBAIJ, 1397 /* 49*/ 0, 1398 MatGetRowIJ_SeqSBAIJ, 1399 MatRestoreRowIJ_SeqSBAIJ, 1400 0, 1401 0, 1402 /* 54*/ 0, 1403 0, 1404 0, 1405 0, 1406 MatSetValuesBlocked_SeqSBAIJ, 1407 /* 59*/ MatGetSubMatrix_SeqSBAIJ, 1408 0, 1409 0, 1410 0, 1411 0, 1412 /* 64*/ 0, 1413 0, 1414 0, 1415 0, 1416 0, 1417 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1418 0, 1419 0, 1420 0, 1421 0, 1422 /* 74*/ 0, 1423 0, 1424 0, 1425 0, 1426 0, 1427 /* 79*/ 0, 1428 0, 1429 0, 1430 MatGetInertia_SeqSBAIJ, 1431 MatLoad_SeqSBAIJ, 1432 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1433 MatIsHermitian_SeqSBAIJ, 1434 MatIsStructurallySymmetric_SeqSBAIJ, 1435 0, 1436 0, 1437 /* 89*/ 0, 1438 0, 1439 0, 1440 0, 1441 0, 1442 /* 94*/ 0, 1443 0, 1444 0, 1445 0, 1446 0, 1447 /* 99*/ 0, 1448 0, 1449 0, 1450 0, 1451 0, 1452 /*104*/ 0, 1453 MatRealPart_SeqSBAIJ, 1454 MatImaginaryPart_SeqSBAIJ, 1455 MatGetRowUpperTriangular_SeqSBAIJ, 1456 MatRestoreRowUpperTriangular_SeqSBAIJ, 1457 /*109*/ 0, 1458 0, 1459 0, 1460 0, 1461 MatMissingDiagonal_SeqSBAIJ, 1462 /*114*/ 0, 1463 0, 1464 0, 1465 0, 1466 0, 1467 /*119*/ 0, 1468 0, 1469 0, 1470 0, 1471 0, 1472 /*124*/ 0, 1473 0, 1474 0, 1475 0, 1476 0, 1477 /*129*/ 0, 1478 0, 1479 0, 1480 0, 1481 0, 1482 /*134*/ 0, 1483 0, 1484 0, 1485 0, 1486 0, 1487 /*139*/ 0, 1488 0, 1489 0, 1490 0, 1491 0, 1492 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ 1493 }; 1494 1495 #undef __FUNCT__ 1496 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1497 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1498 { 1499 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1500 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1501 PetscErrorCode ierr; 1502 1503 PetscFunctionBegin; 1504 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1505 1506 /* allocate space for values if not already there */ 1507 if (!aij->saved_values) { 1508 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 1509 } 1510 1511 /* copy values over */ 1512 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1513 PetscFunctionReturn(0); 1514 } 1515 1516 #undef __FUNCT__ 1517 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1518 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1519 { 1520 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1521 PetscErrorCode ierr; 1522 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1523 1524 PetscFunctionBegin; 1525 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1526 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1527 1528 /* copy values over */ 1529 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1530 PetscFunctionReturn(0); 1531 } 1532 1533 #undef __FUNCT__ 1534 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1535 PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1536 { 1537 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1538 PetscErrorCode ierr; 1539 PetscInt i,mbs,nbs,bs2; 1540 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1541 1542 PetscFunctionBegin; 1543 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1544 B->preallocated = PETSC_TRUE; 1545 1546 ierr = MatSetBlockSize(B,PetscAbs(bs));CHKERRQ(ierr); 1547 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1548 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1549 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1550 1551 mbs = B->rmap->N/bs; 1552 nbs = B->cmap->n/bs; 1553 bs2 = bs*bs; 1554 1555 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"); 1556 1557 if (nz == MAT_SKIP_ALLOCATION) { 1558 skipallocation = PETSC_TRUE; 1559 nz = 0; 1560 } 1561 1562 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1563 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1564 if (nnz) { 1565 for (i=0; i<mbs; i++) { 1566 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]); 1567 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); 1568 } 1569 } 1570 1571 B->ops->mult = MatMult_SeqSBAIJ_N; 1572 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1573 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1574 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1575 1576 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1577 if (!flg) { 1578 switch (bs) { 1579 case 1: 1580 B->ops->mult = MatMult_SeqSBAIJ_1; 1581 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1582 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1583 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1584 break; 1585 case 2: 1586 B->ops->mult = MatMult_SeqSBAIJ_2; 1587 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1588 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1589 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1590 break; 1591 case 3: 1592 B->ops->mult = MatMult_SeqSBAIJ_3; 1593 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1594 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1595 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1596 break; 1597 case 4: 1598 B->ops->mult = MatMult_SeqSBAIJ_4; 1599 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1600 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1601 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1602 break; 1603 case 5: 1604 B->ops->mult = MatMult_SeqSBAIJ_5; 1605 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1606 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1607 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1608 break; 1609 case 6: 1610 B->ops->mult = MatMult_SeqSBAIJ_6; 1611 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1612 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1613 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1614 break; 1615 case 7: 1616 B->ops->mult = MatMult_SeqSBAIJ_7; 1617 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1618 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1619 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1620 break; 1621 } 1622 } 1623 1624 b->mbs = mbs; 1625 b->nbs = nbs; 1626 if (!skipallocation) { 1627 if (!b->imax) { 1628 ierr = PetscMalloc2(mbs,&b->imax,mbs,&b->ilen);CHKERRQ(ierr); 1629 1630 b->free_imax_ilen = PETSC_TRUE; 1631 1632 ierr = PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1633 } 1634 if (!nnz) { 1635 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1636 else if (nz <= 0) nz = 1; 1637 for (i=0; i<mbs; i++) b->imax[i] = nz; 1638 nz = nz*mbs; /* total nz */ 1639 } else { 1640 nz = 0; 1641 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1642 } 1643 /* b->ilen will count nonzeros in each block row so far. */ 1644 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1645 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1646 1647 /* allocate the matrix space */ 1648 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1649 ierr = PetscMalloc3(bs2*nz,&b->a,nz,&b->j,B->rmap->N+1,&b->i);CHKERRQ(ierr); 1650 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1651 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1652 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1653 1654 b->singlemalloc = PETSC_TRUE; 1655 1656 /* pointer to beginning of each row */ 1657 b->i[0] = 0; 1658 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1659 1660 b->free_a = PETSC_TRUE; 1661 b->free_ij = PETSC_TRUE; 1662 } else { 1663 b->free_a = PETSC_FALSE; 1664 b->free_ij = PETSC_FALSE; 1665 } 1666 1667 B->rmap->bs = bs; 1668 b->bs2 = bs2; 1669 b->nz = 0; 1670 b->maxnz = nz; 1671 1672 b->inew = 0; 1673 b->jnew = 0; 1674 b->anew = 0; 1675 b->a2anew = 0; 1676 b->permute = PETSC_FALSE; 1677 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1678 PetscFunctionReturn(0); 1679 } 1680 1681 #undef __FUNCT__ 1682 #define __FUNCT__ "MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ" 1683 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1684 { 1685 PetscInt i,j,m,nz,nz_max=0,*nnz; 1686 PetscScalar *values=0; 1687 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1688 PetscErrorCode ierr; 1689 PetscFunctionBegin; 1690 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs); 1691 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1692 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1693 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1694 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1695 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1696 m = B->rmap->n/bs; 1697 1698 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]); 1699 ierr = PetscMalloc1(m+1,&nnz);CHKERRQ(ierr); 1700 for (i=0; i<m; i++) { 1701 nz = ii[i+1] - ii[i]; 1702 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz); 1703 nz_max = PetscMax(nz_max,nz); 1704 nnz[i] = nz; 1705 } 1706 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1707 ierr = PetscFree(nnz);CHKERRQ(ierr); 1708 1709 values = (PetscScalar*)V; 1710 if (!values) { 1711 ierr = PetscCalloc1(bs*bs*nz_max,&values);CHKERRQ(ierr); 1712 } 1713 for (i=0; i<m; i++) { 1714 PetscInt ncols = ii[i+1] - ii[i]; 1715 const PetscInt *icols = jj + ii[i]; 1716 if (!roworiented || bs == 1) { 1717 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1718 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr); 1719 } else { 1720 for (j=0; j<ncols; j++) { 1721 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1722 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr); 1723 } 1724 } 1725 } 1726 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); } 1727 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1728 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1729 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1730 PetscFunctionReturn(0); 1731 } 1732 1733 /* 1734 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1735 */ 1736 #undef __FUNCT__ 1737 #define __FUNCT__ "MatSeqSBAIJSetNumericFactorization_inplace" 1738 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1739 { 1740 PetscErrorCode ierr; 1741 PetscBool flg = PETSC_FALSE; 1742 PetscInt bs = B->rmap->bs; 1743 1744 PetscFunctionBegin; 1745 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1746 if (flg) bs = 8; 1747 1748 if (!natural) { 1749 switch (bs) { 1750 case 1: 1751 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1752 break; 1753 case 2: 1754 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1755 break; 1756 case 3: 1757 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1758 break; 1759 case 4: 1760 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1761 break; 1762 case 5: 1763 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1764 break; 1765 case 6: 1766 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1767 break; 1768 case 7: 1769 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1770 break; 1771 default: 1772 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1773 break; 1774 } 1775 } else { 1776 switch (bs) { 1777 case 1: 1778 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1779 break; 1780 case 2: 1781 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1782 break; 1783 case 3: 1784 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1785 break; 1786 case 4: 1787 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1788 break; 1789 case 5: 1790 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1791 break; 1792 case 6: 1793 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1794 break; 1795 case 7: 1796 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1797 break; 1798 default: 1799 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1800 break; 1801 } 1802 } 1803 PetscFunctionReturn(0); 1804 } 1805 1806 PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1807 PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1808 1809 #undef __FUNCT__ 1810 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1811 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1812 { 1813 PetscInt n = A->rmap->n; 1814 PetscErrorCode ierr; 1815 1816 PetscFunctionBegin; 1817 #if defined(PETSC_USE_COMPLEX) 1818 if (A->hermitian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported"); 1819 #endif 1820 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 1821 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1822 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1823 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1824 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1825 1826 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1827 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1828 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1829 (*B)->factortype = ftype; 1830 PetscFunctionReturn(0); 1831 } 1832 1833 /*MC 1834 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1835 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1836 1837 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1838 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1839 1840 Options Database Keys: 1841 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1842 1843 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1844 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1845 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. 1846 1847 1848 Level: beginner 1849 1850 .seealso: MatCreateSeqSBAIJ 1851 M*/ 1852 1853 PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1854 1855 #undef __FUNCT__ 1856 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1857 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1858 { 1859 Mat_SeqSBAIJ *b; 1860 PetscErrorCode ierr; 1861 PetscMPIInt size; 1862 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1863 1864 PetscFunctionBegin; 1865 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 1866 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1867 1868 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 1869 B->data = (void*)b; 1870 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1871 1872 B->ops->destroy = MatDestroy_SeqSBAIJ; 1873 B->ops->view = MatView_SeqSBAIJ; 1874 b->row = 0; 1875 b->icol = 0; 1876 b->reallocs = 0; 1877 b->saved_values = 0; 1878 b->inode.limit = 5; 1879 b->inode.max_limit = 5; 1880 1881 b->roworiented = PETSC_TRUE; 1882 b->nonew = 0; 1883 b->diag = 0; 1884 b->solve_work = 0; 1885 b->mult_work = 0; 1886 B->spptr = 0; 1887 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1888 b->keepnonzeropattern = PETSC_FALSE; 1889 1890 b->inew = 0; 1891 b->jnew = 0; 1892 b->anew = 0; 1893 b->a2anew = 0; 1894 b->permute = PETSC_FALSE; 1895 1896 b->ignore_ltriangular = PETSC_TRUE; 1897 1898 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1899 1900 b->getrow_utriangular = PETSC_FALSE; 1901 1902 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 1903 1904 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1905 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1906 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1907 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1908 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1909 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1910 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr); 1911 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqsbstrm_C",MatConvert_SeqSBAIJ_SeqSBSTRM);CHKERRQ(ierr); 1912 #if defined(PETSC_HAVE_ELEMENTAL) 1913 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_elemental_C",MatConvert_SeqSBAIJ_Elemental);CHKERRQ(ierr); 1914 #endif 1915 1916 B->symmetric = PETSC_TRUE; 1917 B->structurally_symmetric = PETSC_TRUE; 1918 B->symmetric_set = PETSC_TRUE; 1919 B->structurally_symmetric_set = PETSC_TRUE; 1920 1921 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1922 1923 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1924 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 1925 if (no_unroll) { 1926 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 1927 } 1928 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 1929 if (no_inode) { 1930 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 1931 } 1932 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 1933 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1934 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1935 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1936 PetscFunctionReturn(0); 1937 } 1938 1939 #undef __FUNCT__ 1940 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1941 /*@C 1942 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1943 compressed row) format. For good matrix assembly performance the 1944 user should preallocate the matrix storage by setting the parameter nz 1945 (or the array nnz). By setting these parameters accurately, performance 1946 during matrix assembly can be increased by more than a factor of 50. 1947 1948 Collective on Mat 1949 1950 Input Parameters: 1951 + B - the symmetric matrix 1952 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 1953 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 1954 . nz - number of block nonzeros per block row (same for all rows) 1955 - nnz - array containing the number of block nonzeros in the upper triangular plus 1956 diagonal portion of each block (possibly different for each block row) or NULL 1957 1958 Options Database Keys: 1959 . -mat_no_unroll - uses code that does not unroll the loops in the 1960 block calculations (much slower) 1961 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1962 1963 Level: intermediate 1964 1965 Notes: 1966 Specify the preallocated storage with either nz or nnz (not both). 1967 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 1968 allocation. See Users-Manual: ch_mat for details. 1969 1970 You can call MatGetInfo() to get information on how effective the preallocation was; 1971 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1972 You can also run with the option -info and look for messages with the string 1973 malloc in them to see if additional memory allocation was needed. 1974 1975 If the nnz parameter is given then the nz parameter is ignored 1976 1977 1978 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 1979 @*/ 1980 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1981 { 1982 PetscErrorCode ierr; 1983 1984 PetscFunctionBegin; 1985 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1986 PetscValidType(B,1); 1987 PetscValidLogicalCollectiveInt(B,bs,2); 1988 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 1989 PetscFunctionReturn(0); 1990 } 1991 1992 #undef __FUNCT__ 1993 #define __FUNCT__ "MatSeqSBAIJSetPreallocationCSR" 1994 /*@C 1995 MatSeqSBAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in symmetric block AIJ format. 1996 1997 Input Parameters: 1998 + B - the matrix 1999 . i - the indices into j for the start of each local row (starts with zero) 2000 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2001 - v - optional values in the matrix 2002 2003 Level: developer 2004 2005 Notes: 2006 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2007 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2008 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2009 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2010 block column and the second index is over columns within a block. 2011 2012 .keywords: matrix, block, aij, compressed row, sparse 2013 2014 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ 2015 @*/ 2016 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2017 { 2018 PetscErrorCode ierr; 2019 2020 PetscFunctionBegin; 2021 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2022 PetscValidType(B,1); 2023 PetscValidLogicalCollectiveInt(B,bs,2); 2024 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2025 PetscFunctionReturn(0); 2026 } 2027 2028 #undef __FUNCT__ 2029 #define __FUNCT__ "MatCreateSeqSBAIJ" 2030 /*@C 2031 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2032 compressed row) format. For good matrix assembly performance the 2033 user should preallocate the matrix storage by setting the parameter nz 2034 (or the array nnz). By setting these parameters accurately, performance 2035 during matrix assembly can be increased by more than a factor of 50. 2036 2037 Collective on MPI_Comm 2038 2039 Input Parameters: 2040 + comm - MPI communicator, set to PETSC_COMM_SELF 2041 . bs - size of block, the blocks are ALWAYS square. One can use MatSetBlockSizes() to set a different row and column blocksize but the row 2042 blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with MatCreateVecs() 2043 . m - number of rows, or number of columns 2044 . nz - number of block nonzeros per block row (same for all rows) 2045 - nnz - array containing the number of block nonzeros in the upper triangular plus 2046 diagonal portion of each block (possibly different for each block row) or NULL 2047 2048 Output Parameter: 2049 . A - the symmetric matrix 2050 2051 Options Database Keys: 2052 . -mat_no_unroll - uses code that does not unroll the loops in the 2053 block calculations (much slower) 2054 . -mat_block_size - size of the blocks to use 2055 2056 Level: intermediate 2057 2058 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2059 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2060 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2061 2062 Notes: 2063 The number of rows and columns must be divisible by blocksize. 2064 This matrix type does not support complex Hermitian operation. 2065 2066 Specify the preallocated storage with either nz or nnz (not both). 2067 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2068 allocation. See Users-Manual: ch_mat for details. 2069 2070 If the nnz parameter is given then the nz parameter is ignored 2071 2072 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2073 @*/ 2074 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2075 { 2076 PetscErrorCode ierr; 2077 2078 PetscFunctionBegin; 2079 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2080 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2081 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2082 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2083 PetscFunctionReturn(0); 2084 } 2085 2086 #undef __FUNCT__ 2087 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2088 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2089 { 2090 Mat C; 2091 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2092 PetscErrorCode ierr; 2093 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2094 2095 PetscFunctionBegin; 2096 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2097 2098 *B = 0; 2099 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2100 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2101 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2102 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2103 c = (Mat_SeqSBAIJ*)C->data; 2104 2105 C->preallocated = PETSC_TRUE; 2106 C->factortype = A->factortype; 2107 c->row = 0; 2108 c->icol = 0; 2109 c->saved_values = 0; 2110 c->keepnonzeropattern = a->keepnonzeropattern; 2111 C->assembled = PETSC_TRUE; 2112 2113 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2114 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2115 c->bs2 = a->bs2; 2116 c->mbs = a->mbs; 2117 c->nbs = a->nbs; 2118 2119 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2120 c->imax = a->imax; 2121 c->ilen = a->ilen; 2122 c->free_imax_ilen = PETSC_FALSE; 2123 } else { 2124 ierr = PetscMalloc2((mbs+1),&c->imax,(mbs+1),&c->ilen);CHKERRQ(ierr); 2125 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2126 for (i=0; i<mbs; i++) { 2127 c->imax[i] = a->imax[i]; 2128 c->ilen[i] = a->ilen[i]; 2129 } 2130 c->free_imax_ilen = PETSC_TRUE; 2131 } 2132 2133 /* allocate the matrix space */ 2134 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2135 ierr = PetscMalloc1(bs2*nz,&c->a);CHKERRQ(ierr); 2136 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2137 c->i = a->i; 2138 c->j = a->j; 2139 c->singlemalloc = PETSC_FALSE; 2140 c->free_a = PETSC_TRUE; 2141 c->free_ij = PETSC_FALSE; 2142 c->parent = A; 2143 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2144 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2145 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2146 } else { 2147 ierr = PetscMalloc3(bs2*nz,&c->a,nz,&c->j,mbs+1,&c->i);CHKERRQ(ierr); 2148 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2149 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2150 c->singlemalloc = PETSC_TRUE; 2151 c->free_a = PETSC_TRUE; 2152 c->free_ij = PETSC_TRUE; 2153 } 2154 if (mbs > 0) { 2155 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2156 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2157 } 2158 if (cpvalues == MAT_COPY_VALUES) { 2159 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2160 } else { 2161 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2162 } 2163 if (a->jshort) { 2164 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2165 /* if the parent matrix is reassembled, this child matrix will never notice */ 2166 ierr = PetscMalloc1(nz,&c->jshort);CHKERRQ(ierr); 2167 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2168 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2169 2170 c->free_jshort = PETSC_TRUE; 2171 } 2172 } 2173 2174 c->roworiented = a->roworiented; 2175 c->nonew = a->nonew; 2176 2177 if (a->diag) { 2178 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2179 c->diag = a->diag; 2180 c->free_diag = PETSC_FALSE; 2181 } else { 2182 ierr = PetscMalloc1(mbs,&c->diag);CHKERRQ(ierr); 2183 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2184 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2185 c->free_diag = PETSC_TRUE; 2186 } 2187 } 2188 c->nz = a->nz; 2189 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2190 c->solve_work = 0; 2191 c->mult_work = 0; 2192 2193 *B = C; 2194 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2195 PetscFunctionReturn(0); 2196 } 2197 2198 #undef __FUNCT__ 2199 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2200 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2201 { 2202 Mat_SeqSBAIJ *a; 2203 PetscErrorCode ierr; 2204 int fd; 2205 PetscMPIInt size; 2206 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs = newmat->rmap->bs; 2207 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2208 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2209 PetscInt *masked,nmask,tmp,bs2,ishift; 2210 PetscScalar *aa; 2211 MPI_Comm comm; 2212 2213 PetscFunctionBegin; 2214 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2215 ierr = PetscOptionsGetInt(((PetscObject)newmat)->prefix,"-matload_block_size",&bs,NULL);CHKERRQ(ierr); 2216 if (bs < 0) bs = 1; 2217 bs2 = bs*bs; 2218 2219 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2220 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2221 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2222 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2223 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2224 M = header[1]; N = header[2]; nz = header[3]; 2225 2226 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2227 2228 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2229 2230 /* 2231 This code adds extra rows to make sure the number of rows is 2232 divisible by the blocksize 2233 */ 2234 mbs = M/bs; 2235 extra_rows = bs - M + bs*(mbs); 2236 if (extra_rows == bs) extra_rows = 0; 2237 else mbs++; 2238 if (extra_rows) { 2239 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2240 } 2241 2242 /* Set global sizes if not already set */ 2243 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2244 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2245 } else { /* Check if the matrix global sizes are correct */ 2246 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2247 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); 2248 } 2249 2250 /* read in row lengths */ 2251 ierr = PetscMalloc1(M+extra_rows,&rowlengths);CHKERRQ(ierr); 2252 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2253 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2254 2255 /* read in column indices */ 2256 ierr = PetscMalloc1(nz+extra_rows,&jj);CHKERRQ(ierr); 2257 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2258 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2259 2260 /* loop over row lengths determining block row lengths */ 2261 ierr = PetscCalloc1(mbs,&s_browlengths);CHKERRQ(ierr); 2262 ierr = PetscMalloc2(mbs,&mask,mbs,&masked);CHKERRQ(ierr); 2263 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2264 rowcount = 0; 2265 nzcount = 0; 2266 for (i=0; i<mbs; i++) { 2267 nmask = 0; 2268 for (j=0; j<bs; j++) { 2269 kmax = rowlengths[rowcount]; 2270 for (k=0; k<kmax; k++) { 2271 tmp = jj[nzcount++]/bs; /* block col. index */ 2272 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2273 } 2274 rowcount++; 2275 } 2276 s_browlengths[i] += nmask; 2277 2278 /* zero out the mask elements we set */ 2279 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2280 } 2281 2282 /* Do preallocation */ 2283 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2284 a = (Mat_SeqSBAIJ*)newmat->data; 2285 2286 /* set matrix "i" values */ 2287 a->i[0] = 0; 2288 for (i=1; i<= mbs; i++) { 2289 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2290 a->ilen[i-1] = s_browlengths[i-1]; 2291 } 2292 a->nz = a->i[mbs]; 2293 2294 /* read in nonzero values */ 2295 ierr = PetscMalloc1(nz+extra_rows,&aa);CHKERRQ(ierr); 2296 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2297 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2298 2299 /* set "a" and "j" values into matrix */ 2300 nzcount = 0; jcount = 0; 2301 for (i=0; i<mbs; i++) { 2302 nzcountb = nzcount; 2303 nmask = 0; 2304 for (j=0; j<bs; j++) { 2305 kmax = rowlengths[i*bs+j]; 2306 for (k=0; k<kmax; k++) { 2307 tmp = jj[nzcount++]/bs; /* block col. index */ 2308 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2309 } 2310 } 2311 /* sort the masked values */ 2312 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2313 2314 /* set "j" values into matrix */ 2315 maskcount = 1; 2316 for (j=0; j<nmask; j++) { 2317 a->j[jcount++] = masked[j]; 2318 mask[masked[j]] = maskcount++; 2319 } 2320 2321 /* set "a" values into matrix */ 2322 ishift = bs2*a->i[i]; 2323 for (j=0; j<bs; j++) { 2324 kmax = rowlengths[i*bs+j]; 2325 for (k=0; k<kmax; k++) { 2326 tmp = jj[nzcountb]/bs; /* block col. index */ 2327 if (tmp >= i) { 2328 block = mask[tmp] - 1; 2329 point = jj[nzcountb] - bs*tmp; 2330 idx = ishift + bs2*block + j + bs*point; 2331 a->a[idx] = aa[nzcountb]; 2332 } 2333 nzcountb++; 2334 } 2335 } 2336 /* zero out the mask elements we set */ 2337 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2338 } 2339 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2340 2341 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2342 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2343 ierr = PetscFree(aa);CHKERRQ(ierr); 2344 ierr = PetscFree(jj);CHKERRQ(ierr); 2345 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2346 2347 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2348 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2349 PetscFunctionReturn(0); 2350 } 2351 2352 #undef __FUNCT__ 2353 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2354 /*@ 2355 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2356 (upper triangular entries in CSR format) provided by the user. 2357 2358 Collective on MPI_Comm 2359 2360 Input Parameters: 2361 + comm - must be an MPI communicator of size 1 2362 . bs - size of block 2363 . m - number of rows 2364 . n - number of columns 2365 . i - row indices 2366 . j - column indices 2367 - a - matrix values 2368 2369 Output Parameter: 2370 . mat - the matrix 2371 2372 Level: advanced 2373 2374 Notes: 2375 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2376 once the matrix is destroyed 2377 2378 You cannot set new nonzero locations into this matrix, that will generate an error. 2379 2380 The i and j indices are 0 based 2381 2382 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 2383 it is the regular CSR format excluding the lower triangular elements. 2384 2385 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2386 2387 @*/ 2388 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 2389 { 2390 PetscErrorCode ierr; 2391 PetscInt ii; 2392 Mat_SeqSBAIJ *sbaij; 2393 2394 PetscFunctionBegin; 2395 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2396 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2397 2398 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2399 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2400 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2401 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2402 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2403 ierr = PetscMalloc2(m,&sbaij->imax,m,&sbaij->ilen);CHKERRQ(ierr); 2404 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2405 2406 sbaij->i = i; 2407 sbaij->j = j; 2408 sbaij->a = a; 2409 2410 sbaij->singlemalloc = PETSC_FALSE; 2411 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2412 sbaij->free_a = PETSC_FALSE; 2413 sbaij->free_ij = PETSC_FALSE; 2414 2415 for (ii=0; ii<m; ii++) { 2416 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2417 #if defined(PETSC_USE_DEBUG) 2418 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]); 2419 #endif 2420 } 2421 #if defined(PETSC_USE_DEBUG) 2422 for (ii=0; ii<sbaij->i[m]; ii++) { 2423 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2424 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]); 2425 } 2426 #endif 2427 2428 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2429 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2430 PetscFunctionReturn(0); 2431 } 2432 2433 #undef __FUNCT__ 2434 #define __FUNCT__ "MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ" 2435 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqSBAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2436 { 2437 PetscErrorCode ierr; 2438 2439 PetscFunctionBegin; 2440 ierr = MatCreateMPIMatConcatenateSeqMat_MPISBAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2441 PetscFunctionReturn(0); 2442 } 2443 2444 2445 2446 2447