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