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