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