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