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