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 1117 PetscFunctionBegin; 1118 nz = baij->maxnz; 1119 n = mat->cmap->n; 1120 for (i=0; i<nz; i++) { 1121 baij->j[i] = indices[i]; 1122 } 1123 baij->nz = nz; 1124 for (i=0; i<n; i++) { 1125 baij->ilen[i] = baij->imax[i]; 1126 } 1127 PetscFunctionReturn(0); 1128 } 1129 EXTERN_C_END 1130 1131 #undef __FUNCT__ 1132 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1133 /*@ 1134 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1135 in the matrix. 1136 1137 Input Parameters: 1138 + mat - the SeqSBAIJ matrix 1139 - indices - the column indices 1140 1141 Level: advanced 1142 1143 Notes: 1144 This can be called if you have precomputed the nonzero structure of the 1145 matrix and want to provide it to the matrix object to improve the performance 1146 of the MatSetValues() operation. 1147 1148 You MUST have set the correct numbers of nonzeros per row in the call to 1149 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1150 1151 MUST be called before any calls to MatSetValues() 1152 1153 .seealso: MatCreateSeqSBAIJ 1154 @*/ 1155 PetscErrorCode MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1156 { 1157 PetscErrorCode ierr; 1158 1159 PetscFunctionBegin; 1160 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1161 PetscValidPointer(indices,2); 1162 ierr = PetscUseMethod(mat,"MatSeqSBAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1168 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1169 { 1170 PetscErrorCode ierr; 1171 1172 PetscFunctionBegin; 1173 /* If the two matrices have the same copy implementation, use fast copy. */ 1174 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1175 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1176 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1177 1178 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"); 1179 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1180 } else { 1181 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1182 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1183 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1184 } 1185 PetscFunctionReturn(0); 1186 } 1187 1188 #undef __FUNCT__ 1189 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1190 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1191 { 1192 PetscErrorCode ierr; 1193 1194 PetscFunctionBegin; 1195 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,-PetscMax(A->rmap->bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1196 PetscFunctionReturn(0); 1197 } 1198 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1201 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1202 { 1203 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1204 PetscFunctionBegin; 1205 *array = a->a; 1206 PetscFunctionReturn(0); 1207 } 1208 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1211 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1212 { 1213 PetscFunctionBegin; 1214 PetscFunctionReturn(0); 1215 } 1216 1217 #undef __FUNCT__ 1218 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1219 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1220 { 1221 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1222 PetscErrorCode ierr; 1223 PetscInt i,bs=Y->rmap->bs,bs2=bs*bs,j; 1224 PetscBLASInt one = 1; 1225 1226 PetscFunctionBegin; 1227 if (str == SAME_NONZERO_PATTERN) { 1228 PetscScalar alpha = a; 1229 PetscBLASInt bnz = PetscBLASIntCast(x->nz*bs2); 1230 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1231 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1232 if (y->xtoy && y->XtoY != X) { 1233 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1234 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 1235 } 1236 if (!y->xtoy) { /* get xtoy */ 1237 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1238 y->XtoY = X; 1239 } 1240 for (i=0; i<x->nz; i++) { 1241 j = 0; 1242 while (j < bs2){ 1243 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1244 j++; 1245 } 1246 } 1247 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); 1248 } else { 1249 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1250 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1251 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1252 } 1253 PetscFunctionReturn(0); 1254 } 1255 1256 #undef __FUNCT__ 1257 #define __FUNCT__ "MatSetBlockSize_SeqSBAIJ" 1258 PetscErrorCode MatSetBlockSize_SeqSBAIJ(Mat A,PetscInt bs) 1259 { 1260 PetscInt rbs,cbs; 1261 PetscErrorCode ierr; 1262 1263 PetscFunctionBegin; 1264 ierr = PetscLayoutGetBlockSize(A->rmap,&rbs);CHKERRQ(ierr); 1265 ierr = PetscLayoutGetBlockSize(A->cmap,&cbs);CHKERRQ(ierr); 1266 if (rbs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Attempt to set block size %d with SBAIJ %d",bs,rbs); 1267 if (cbs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Attempt to set block size %d with SBAIJ %d",bs,cbs); 1268 PetscFunctionReturn(0); 1269 } 1270 1271 #undef __FUNCT__ 1272 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1273 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1274 { 1275 PetscFunctionBegin; 1276 *flg = PETSC_TRUE; 1277 PetscFunctionReturn(0); 1278 } 1279 1280 #undef __FUNCT__ 1281 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1282 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscBool *flg) 1283 { 1284 PetscFunctionBegin; 1285 *flg = PETSC_TRUE; 1286 PetscFunctionReturn(0); 1287 } 1288 1289 #undef __FUNCT__ 1290 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1291 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscBool *flg) 1292 { 1293 PetscFunctionBegin; 1294 *flg = PETSC_FALSE; 1295 PetscFunctionReturn(0); 1296 } 1297 1298 #undef __FUNCT__ 1299 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1300 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1301 { 1302 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1303 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1304 MatScalar *aa = a->a; 1305 1306 PetscFunctionBegin; 1307 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1308 PetscFunctionReturn(0); 1309 } 1310 1311 #undef __FUNCT__ 1312 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1313 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1314 { 1315 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1316 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1317 MatScalar *aa = a->a; 1318 1319 PetscFunctionBegin; 1320 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1321 PetscFunctionReturn(0); 1322 } 1323 1324 #undef __FUNCT__ 1325 #define __FUNCT__ "MatZeroRowsColumns_SeqSBAIJ" 1326 PetscErrorCode MatZeroRowsColumns_SeqSBAIJ(Mat A,PetscInt is_n,const PetscInt is_idx[],PetscScalar diag,Vec x, Vec b) 1327 { 1328 Mat_SeqSBAIJ *baij=(Mat_SeqSBAIJ*)A->data; 1329 PetscErrorCode ierr; 1330 PetscInt i,j,k,count; 1331 PetscInt bs=A->rmap->bs,bs2=baij->bs2,row,col; 1332 PetscScalar zero = 0.0; 1333 MatScalar *aa; 1334 const PetscScalar *xx; 1335 PetscScalar *bb; 1336 PetscBool *zeroed,vecs = PETSC_FALSE; 1337 1338 PetscFunctionBegin; 1339 /* fix right hand side if needed */ 1340 if (x && b) { 1341 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1342 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1343 vecs = PETSC_TRUE; 1344 } 1345 A->same_nonzero = PETSC_TRUE; 1346 1347 /* zero the columns */ 1348 ierr = PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);CHKERRQ(ierr); 1349 ierr = PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));CHKERRQ(ierr); 1350 for (i=0; i<is_n; i++) { 1351 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]); 1352 zeroed[is_idx[i]] = PETSC_TRUE; 1353 } 1354 if (vecs) { 1355 for (i=0; i<A->rmap->N; i++) { 1356 row = i/bs; 1357 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1358 for (k=0; k<bs; k++) { 1359 col = bs*baij->j[j] + k; 1360 if (col <= i) continue; 1361 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1362 if (!zeroed[i] && zeroed[col]) { 1363 bb[i] -= aa[0]*xx[col]; 1364 } 1365 if (zeroed[i] && !zeroed[col]) { 1366 bb[col] -= aa[0]*xx[i]; 1367 } 1368 } 1369 } 1370 } 1371 for (i=0; i<is_n; i++) { 1372 bb[is_idx[i]] = diag*xx[is_idx[i]]; 1373 } 1374 } 1375 1376 for (i=0; i<A->rmap->N; i++) { 1377 if (!zeroed[i]) { 1378 row = i/bs; 1379 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1380 for (k=0; k<bs; k++) { 1381 col = bs*baij->j[j] + k; 1382 if (zeroed[col]) { 1383 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1384 aa[0] = 0.0; 1385 } 1386 } 1387 } 1388 } 1389 } 1390 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1391 if (vecs) { 1392 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1393 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1394 } 1395 1396 /* zero the rows */ 1397 for (i=0; i<is_n; i++) { 1398 row = is_idx[i]; 1399 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1400 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1401 for (k=0; k<count; k++) { 1402 aa[0] = zero; 1403 aa += bs; 1404 } 1405 if (diag != 0.0) { 1406 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1407 } 1408 } 1409 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1410 PetscFunctionReturn(0); 1411 } 1412 1413 /* -------------------------------------------------------------------*/ 1414 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1415 MatGetRow_SeqSBAIJ, 1416 MatRestoreRow_SeqSBAIJ, 1417 MatMult_SeqSBAIJ_N, 1418 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1419 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1420 MatMultAdd_SeqSBAIJ_N, 1421 0, 1422 0, 1423 0, 1424 /*10*/ 0, 1425 0, 1426 MatCholeskyFactor_SeqSBAIJ, 1427 MatSOR_SeqSBAIJ, 1428 MatTranspose_SeqSBAIJ, 1429 /*15*/ MatGetInfo_SeqSBAIJ, 1430 MatEqual_SeqSBAIJ, 1431 MatGetDiagonal_SeqSBAIJ, 1432 MatDiagonalScale_SeqSBAIJ, 1433 MatNorm_SeqSBAIJ, 1434 /*20*/ 0, 1435 MatAssemblyEnd_SeqSBAIJ, 1436 MatSetOption_SeqSBAIJ, 1437 MatZeroEntries_SeqSBAIJ, 1438 /*24*/ 0, 1439 0, 1440 0, 1441 0, 1442 0, 1443 /*29*/ MatSetUpPreallocation_SeqSBAIJ, 1444 0, 1445 0, 1446 MatGetArray_SeqSBAIJ, 1447 MatRestoreArray_SeqSBAIJ, 1448 /*34*/ MatDuplicate_SeqSBAIJ, 1449 0, 1450 0, 1451 0, 1452 MatICCFactor_SeqSBAIJ, 1453 /*39*/ MatAXPY_SeqSBAIJ, 1454 MatGetSubMatrices_SeqSBAIJ, 1455 MatIncreaseOverlap_SeqSBAIJ, 1456 MatGetValues_SeqSBAIJ, 1457 MatCopy_SeqSBAIJ, 1458 /*44*/ 0, 1459 MatScale_SeqSBAIJ, 1460 0, 1461 0, 1462 MatZeroRowsColumns_SeqSBAIJ, 1463 /*49*/ MatSetBlockSize_SeqSBAIJ, 1464 MatGetRowIJ_SeqSBAIJ, 1465 MatRestoreRowIJ_SeqSBAIJ, 1466 0, 1467 0, 1468 /*54*/ 0, 1469 0, 1470 0, 1471 0, 1472 MatSetValuesBlocked_SeqSBAIJ, 1473 /*59*/ MatGetSubMatrix_SeqSBAIJ, 1474 0, 1475 0, 1476 0, 1477 0, 1478 /*64*/ 0, 1479 0, 1480 0, 1481 0, 1482 0, 1483 /*69*/ MatGetRowMaxAbs_SeqSBAIJ, 1484 0, 1485 0, 1486 0, 1487 0, 1488 /*74*/ 0, 1489 0, 1490 0, 1491 0, 1492 0, 1493 /*79*/ 0, 1494 0, 1495 0, 1496 MatGetInertia_SeqSBAIJ, 1497 MatLoad_SeqSBAIJ, 1498 /*84*/ MatIsSymmetric_SeqSBAIJ, 1499 MatIsHermitian_SeqSBAIJ, 1500 MatIsStructurallySymmetric_SeqSBAIJ, 1501 0, 1502 0, 1503 /*89*/ 0, 1504 0, 1505 0, 1506 0, 1507 0, 1508 /*94*/ 0, 1509 0, 1510 0, 1511 0, 1512 0, 1513 /*99*/ 0, 1514 0, 1515 0, 1516 0, 1517 0, 1518 /*104*/0, 1519 MatRealPart_SeqSBAIJ, 1520 MatImaginaryPart_SeqSBAIJ, 1521 MatGetRowUpperTriangular_SeqSBAIJ, 1522 MatRestoreRowUpperTriangular_SeqSBAIJ, 1523 /*109*/0, 1524 0, 1525 0, 1526 0, 1527 MatMissingDiagonal_SeqSBAIJ, 1528 /*114*/0, 1529 0, 1530 0, 1531 0, 1532 0, 1533 /*119*/0, 1534 0, 1535 0, 1536 0 1537 }; 1538 1539 EXTERN_C_BEGIN 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1542 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1543 { 1544 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1545 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1546 PetscErrorCode ierr; 1547 1548 PetscFunctionBegin; 1549 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1550 1551 /* allocate space for values if not already there */ 1552 if (!aij->saved_values) { 1553 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1554 } 1555 1556 /* copy values over */ 1557 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1558 PetscFunctionReturn(0); 1559 } 1560 EXTERN_C_END 1561 1562 EXTERN_C_BEGIN 1563 #undef __FUNCT__ 1564 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1565 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1566 { 1567 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1568 PetscErrorCode ierr; 1569 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1570 1571 PetscFunctionBegin; 1572 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1573 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1574 1575 /* copy values over */ 1576 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1577 PetscFunctionReturn(0); 1578 } 1579 EXTERN_C_END 1580 1581 EXTERN_C_BEGIN 1582 #undef __FUNCT__ 1583 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1584 PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1585 { 1586 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1587 PetscErrorCode ierr; 1588 PetscInt i,mbs,bs2, newbs = PetscAbs(bs); 1589 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE; 1590 1591 PetscFunctionBegin; 1592 B->preallocated = PETSC_TRUE; 1593 if (bs < 0) { 1594 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for MPISBAIJ matrix","Mat");CHKERRQ(ierr); 1595 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSeqSBAIJSetPreallocation",newbs,&newbs,PETSC_NULL);CHKERRQ(ierr); 1596 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1597 bs = PetscAbs(bs); 1598 } 1599 if (nnz && newbs != bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot change blocksize from command line if setting nnz"); 1600 bs = newbs; 1601 1602 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1603 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1604 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1605 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1606 1607 mbs = B->rmap->N/bs; 1608 bs2 = bs*bs; 1609 1610 if (mbs*bs != B->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1611 1612 if (nz == MAT_SKIP_ALLOCATION) { 1613 skipallocation = PETSC_TRUE; 1614 nz = 0; 1615 } 1616 1617 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1618 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1619 if (nnz) { 1620 for (i=0; i<mbs; i++) { 1621 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]); 1622 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); 1623 } 1624 } 1625 1626 B->ops->mult = MatMult_SeqSBAIJ_N; 1627 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1628 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1629 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1630 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1631 if (!flg) { 1632 switch (bs) { 1633 case 1: 1634 B->ops->mult = MatMult_SeqSBAIJ_1; 1635 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1636 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1637 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1638 break; 1639 case 2: 1640 B->ops->mult = MatMult_SeqSBAIJ_2; 1641 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1642 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1643 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1644 break; 1645 case 3: 1646 B->ops->mult = MatMult_SeqSBAIJ_3; 1647 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1648 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1649 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1650 break; 1651 case 4: 1652 B->ops->mult = MatMult_SeqSBAIJ_4; 1653 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1654 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1655 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1656 break; 1657 case 5: 1658 B->ops->mult = MatMult_SeqSBAIJ_5; 1659 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1660 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1661 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1662 break; 1663 case 6: 1664 B->ops->mult = MatMult_SeqSBAIJ_6; 1665 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1666 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1667 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1668 break; 1669 case 7: 1670 B->ops->mult = MatMult_SeqSBAIJ_7; 1671 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1672 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1673 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1674 break; 1675 } 1676 } 1677 1678 b->mbs = mbs; 1679 b->nbs = mbs; 1680 if (!skipallocation) { 1681 if (!b->imax) { 1682 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1683 b->free_imax_ilen = PETSC_TRUE; 1684 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1685 } 1686 if (!nnz) { 1687 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1688 else if (nz <= 0) nz = 1; 1689 for (i=0; i<mbs; i++) { 1690 b->imax[i] = nz; 1691 } 1692 nz = nz*mbs; /* total nz */ 1693 } else { 1694 nz = 0; 1695 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1696 } 1697 /* b->ilen will count nonzeros in each block row so far. */ 1698 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1699 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1700 1701 /* allocate the matrix space */ 1702 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1703 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1704 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1705 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1706 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1707 b->singlemalloc = PETSC_TRUE; 1708 1709 /* pointer to beginning of each row */ 1710 b->i[0] = 0; 1711 for (i=1; i<mbs+1; i++) { 1712 b->i[i] = b->i[i-1] + b->imax[i-1]; 1713 } 1714 b->free_a = PETSC_TRUE; 1715 b->free_ij = PETSC_TRUE; 1716 } else { 1717 b->free_a = PETSC_FALSE; 1718 b->free_ij = PETSC_FALSE; 1719 } 1720 1721 B->rmap->bs = bs; 1722 b->bs2 = bs2; 1723 b->nz = 0; 1724 b->maxnz = nz; 1725 1726 b->inew = 0; 1727 b->jnew = 0; 1728 b->anew = 0; 1729 b->a2anew = 0; 1730 b->permute = PETSC_FALSE; 1731 PetscFunctionReturn(0); 1732 } 1733 EXTERN_C_END 1734 1735 /* 1736 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1737 */ 1738 #undef __FUNCT__ 1739 #define __FUNCT__ "MatSeqSBAIJSetNumericFactorization_inplace" 1740 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1741 { 1742 PetscErrorCode ierr; 1743 PetscBool flg = PETSC_FALSE; 1744 PetscInt bs = B->rmap->bs; 1745 1746 PetscFunctionBegin; 1747 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1748 if (flg) bs = 8; 1749 1750 if (!natural) { 1751 switch (bs) { 1752 case 1: 1753 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1754 break; 1755 case 2: 1756 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1757 break; 1758 case 3: 1759 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1760 break; 1761 case 4: 1762 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1763 break; 1764 case 5: 1765 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1766 break; 1767 case 6: 1768 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1769 break; 1770 case 7: 1771 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1772 break; 1773 default: 1774 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1775 break; 1776 } 1777 } else { 1778 switch (bs) { 1779 case 1: 1780 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1781 break; 1782 case 2: 1783 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1784 break; 1785 case 3: 1786 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1787 break; 1788 case 4: 1789 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1790 break; 1791 case 5: 1792 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1793 break; 1794 case 6: 1795 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1796 break; 1797 case 7: 1798 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1799 break; 1800 default: 1801 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1802 break; 1803 } 1804 } 1805 PetscFunctionReturn(0); 1806 } 1807 1808 EXTERN_C_BEGIN 1809 extern PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1810 extern PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1811 EXTERN_C_END 1812 1813 1814 EXTERN_C_BEGIN 1815 #undef __FUNCT__ 1816 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1817 PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1818 { 1819 PetscInt n = A->rmap->n; 1820 PetscErrorCode ierr; 1821 1822 PetscFunctionBegin; 1823 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1824 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1825 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1826 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1827 ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1828 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1829 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1830 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1831 (*B)->factortype = ftype; 1832 PetscFunctionReturn(0); 1833 } 1834 EXTERN_C_END 1835 1836 EXTERN_C_BEGIN 1837 #undef __FUNCT__ 1838 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1839 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscBool *flg) 1840 { 1841 PetscFunctionBegin; 1842 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1843 *flg = PETSC_TRUE; 1844 } else { 1845 *flg = PETSC_FALSE; 1846 } 1847 PetscFunctionReturn(0); 1848 } 1849 EXTERN_C_END 1850 1851 EXTERN_C_BEGIN 1852 #if defined(PETSC_HAVE_MUMPS) 1853 extern PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 1854 #endif 1855 #if defined(PETSC_HAVE_SPOOLES) 1856 extern PetscErrorCode MatGetFactor_seqsbaij_spooles(Mat,MatFactorType,Mat*); 1857 #endif 1858 #if defined(PETSC_HAVE_PASTIX) 1859 extern PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1860 #endif 1861 #if defined(PETSC_HAVE_CHOLMOD) 1862 extern PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*); 1863 #endif 1864 extern PetscErrorCode MatGetFactor_seqsbaij_sbstrm(Mat,MatFactorType,Mat*); 1865 EXTERN_C_END 1866 1867 /*MC 1868 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1869 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1870 1871 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1872 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1873 1874 Options Database Keys: 1875 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1876 1877 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1878 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1879 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. 1880 1881 1882 Level: beginner 1883 1884 .seealso: MatCreateSeqSBAIJ 1885 M*/ 1886 1887 EXTERN_C_BEGIN 1888 extern PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1889 EXTERN_C_END 1890 1891 1892 EXTERN_C_BEGIN 1893 #undef __FUNCT__ 1894 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1895 PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1896 { 1897 Mat_SeqSBAIJ *b; 1898 PetscErrorCode ierr; 1899 PetscMPIInt size; 1900 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1901 1902 PetscFunctionBegin; 1903 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1904 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1905 1906 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1907 B->data = (void*)b; 1908 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1909 B->ops->destroy = MatDestroy_SeqSBAIJ; 1910 B->ops->view = MatView_SeqSBAIJ; 1911 b->row = 0; 1912 b->icol = 0; 1913 b->reallocs = 0; 1914 b->saved_values = 0; 1915 b->inode.limit = 5; 1916 b->inode.max_limit = 5; 1917 1918 b->roworiented = PETSC_TRUE; 1919 b->nonew = 0; 1920 b->diag = 0; 1921 b->solve_work = 0; 1922 b->mult_work = 0; 1923 B->spptr = 0; 1924 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1925 b->keepnonzeropattern = PETSC_FALSE; 1926 b->xtoy = 0; 1927 b->XtoY = 0; 1928 1929 b->inew = 0; 1930 b->jnew = 0; 1931 b->anew = 0; 1932 b->a2anew = 0; 1933 b->permute = PETSC_FALSE; 1934 1935 b->ignore_ltriangular = PETSC_TRUE; 1936 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,PETSC_NULL);CHKERRQ(ierr); 1937 1938 b->getrow_utriangular = PETSC_FALSE; 1939 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,PETSC_NULL);CHKERRQ(ierr); 1940 1941 #if defined(PETSC_HAVE_PASTIX) 1942 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 1943 "MatGetFactor_seqsbaij_pastix", 1944 MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1945 #endif 1946 #if defined(PETSC_HAVE_SPOOLES) 1947 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 1948 "MatGetFactor_seqsbaij_spooles", 1949 MatGetFactor_seqsbaij_spooles);CHKERRQ(ierr); 1950 #endif 1951 #if defined(PETSC_HAVE_MUMPS) 1952 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 1953 "MatGetFactor_sbaij_mumps", 1954 MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 1955 #endif 1956 #if defined(PETSC_HAVE_CHOLMOD) 1957 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 1958 "MatGetFactor_seqsbaij_cholmod", 1959 MatGetFactor_seqsbaij_cholmod);CHKERRQ(ierr); 1960 #endif 1961 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 1962 "MatGetFactorAvailable_seqsbaij_petsc", 1963 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1964 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 1965 "MatGetFactor_seqsbaij_petsc", 1966 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1967 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_sbstrm_C", 1968 "MatGetFactor_seqsbaij_sbstrm", 1969 MatGetFactor_seqsbaij_sbstrm);CHKERRQ(ierr); 1970 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1971 "MatStoreValues_SeqSBAIJ", 1972 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1973 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1974 "MatRetrieveValues_SeqSBAIJ", 1975 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1976 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1977 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1978 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1979 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1980 "MatConvert_SeqSBAIJ_SeqAIJ", 1981 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1982 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1983 "MatConvert_SeqSBAIJ_SeqBAIJ", 1984 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1985 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1986 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1987 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1988 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqsbstrm_C", 1989 "MatConvert_SeqSBAIJ_SeqSBSTRM", 1990 MatConvert_SeqSBAIJ_SeqSBSTRM);CHKERRQ(ierr); 1991 1992 B->symmetric = PETSC_TRUE; 1993 B->structurally_symmetric = PETSC_TRUE; 1994 B->symmetric_set = PETSC_TRUE; 1995 B->structurally_symmetric_set = PETSC_TRUE; 1996 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1997 1998 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1999 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",PETSC_NULL,no_unroll,&no_unroll,PETSC_NULL);CHKERRQ(ierr); 2000 if (no_unroll) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr);} 2001 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",PETSC_NULL,no_inode,&no_inode,PETSC_NULL);CHKERRQ(ierr); 2002 if (no_inode) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr);} 2003 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); 2004 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2005 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 2006 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 2007 2008 PetscFunctionReturn(0); 2009 } 2010 EXTERN_C_END 2011 2012 #undef __FUNCT__ 2013 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 2014 /*@C 2015 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2016 compressed row) format. For good matrix assembly performance the 2017 user should preallocate the matrix storage by setting the parameter nz 2018 (or the array nnz). By setting these parameters accurately, performance 2019 during matrix assembly can be increased by more than a factor of 50. 2020 2021 Collective on Mat 2022 2023 Input Parameters: 2024 + A - the symmetric matrix 2025 . bs - size of block 2026 . nz - number of block nonzeros per block row (same for all rows) 2027 - nnz - array containing the number of block nonzeros in the upper triangular plus 2028 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 2029 2030 Options Database Keys: 2031 . -mat_no_unroll - uses code that does not unroll the loops in the 2032 block calculations (much slower) 2033 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2034 2035 Level: intermediate 2036 2037 Notes: 2038 Specify the preallocated storage with either nz or nnz (not both). 2039 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2040 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2041 2042 You can call MatGetInfo() to get information on how effective the preallocation was; 2043 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2044 You can also run with the option -info and look for messages with the string 2045 malloc in them to see if additional memory allocation was needed. 2046 2047 If the nnz parameter is given then the nz parameter is ignored 2048 2049 2050 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 2051 @*/ 2052 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2053 { 2054 PetscErrorCode ierr; 2055 2056 PetscFunctionBegin; 2057 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2058 PetscFunctionReturn(0); 2059 } 2060 2061 #undef __FUNCT__ 2062 #define __FUNCT__ "MatCreateSeqSBAIJ" 2063 /*@C 2064 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2065 compressed row) format. For good matrix assembly performance the 2066 user should preallocate the matrix storage by setting the parameter nz 2067 (or the array nnz). By setting these parameters accurately, performance 2068 during matrix assembly can be increased by more than a factor of 50. 2069 2070 Collective on MPI_Comm 2071 2072 Input Parameters: 2073 + comm - MPI communicator, set to PETSC_COMM_SELF 2074 . bs - size of block 2075 . m - number of rows, or number of columns 2076 . nz - number of block nonzeros per block row (same for all rows) 2077 - nnz - array containing the number of block nonzeros in the upper triangular plus 2078 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 2079 2080 Output Parameter: 2081 . A - the symmetric matrix 2082 2083 Options Database Keys: 2084 . -mat_no_unroll - uses code that does not unroll the loops in the 2085 block calculations (much slower) 2086 . -mat_block_size - size of the blocks to use 2087 2088 Level: intermediate 2089 2090 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2091 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2092 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2093 2094 Notes: 2095 The number of rows and columns must be divisible by blocksize. 2096 This matrix type does not support complex Hermitian operation. 2097 2098 Specify the preallocated storage with either nz or nnz (not both). 2099 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2100 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2101 2102 If the nnz parameter is given then the nz parameter is ignored 2103 2104 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 2105 @*/ 2106 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2107 { 2108 PetscErrorCode ierr; 2109 2110 PetscFunctionBegin; 2111 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2112 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2113 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2114 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2115 PetscFunctionReturn(0); 2116 } 2117 2118 #undef __FUNCT__ 2119 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2120 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2121 { 2122 Mat C; 2123 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2124 PetscErrorCode ierr; 2125 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2126 2127 PetscFunctionBegin; 2128 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2129 2130 *B = 0; 2131 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 2132 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2133 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2134 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2135 c = (Mat_SeqSBAIJ*)C->data; 2136 2137 C->preallocated = PETSC_TRUE; 2138 C->factortype = A->factortype; 2139 c->row = 0; 2140 c->icol = 0; 2141 c->saved_values = 0; 2142 c->keepnonzeropattern = a->keepnonzeropattern; 2143 C->assembled = PETSC_TRUE; 2144 2145 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2146 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2147 c->bs2 = a->bs2; 2148 c->mbs = a->mbs; 2149 c->nbs = a->nbs; 2150 2151 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2152 c->imax = a->imax; 2153 c->ilen = a->ilen; 2154 c->free_imax_ilen = PETSC_FALSE; 2155 } else { 2156 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 2157 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2158 for (i=0; i<mbs; i++) { 2159 c->imax[i] = a->imax[i]; 2160 c->ilen[i] = a->ilen[i]; 2161 } 2162 c->free_imax_ilen = PETSC_TRUE; 2163 } 2164 2165 /* allocate the matrix space */ 2166 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2167 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2168 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2169 c->singlemalloc = PETSC_FALSE; 2170 c->free_ij = PETSC_FALSE; 2171 c->parent = A; 2172 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2173 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2174 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2175 } else { 2176 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2177 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2178 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2179 c->singlemalloc = PETSC_TRUE; 2180 c->free_ij = PETSC_TRUE; 2181 } 2182 if (mbs > 0) { 2183 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2184 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2185 } 2186 if (cpvalues == MAT_COPY_VALUES) { 2187 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2188 } else { 2189 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2190 } 2191 if (a->jshort) { 2192 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2193 c->jshort = a->jshort; 2194 c->free_jshort = PETSC_FALSE; 2195 } else { 2196 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2197 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2198 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2199 c->free_jshort = PETSC_TRUE; 2200 } 2201 } 2202 } 2203 2204 c->roworiented = a->roworiented; 2205 c->nonew = a->nonew; 2206 2207 if (a->diag) { 2208 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2209 c->diag = a->diag; 2210 c->free_diag = PETSC_FALSE; 2211 } else { 2212 ierr = PetscMalloc(mbs*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2213 ierr = PetscLogObjectMemory(C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2214 for (i=0; i<mbs; i++) { 2215 c->diag[i] = a->diag[i]; 2216 } 2217 c->free_diag = PETSC_TRUE; 2218 } 2219 } else c->diag = 0; 2220 c->nz = a->nz; 2221 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2222 c->solve_work = 0; 2223 c->mult_work = 0; 2224 c->free_a = PETSC_TRUE; 2225 *B = C; 2226 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2227 PetscFunctionReturn(0); 2228 } 2229 2230 #undef __FUNCT__ 2231 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2232 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2233 { 2234 Mat_SeqSBAIJ *a; 2235 PetscErrorCode ierr; 2236 int fd; 2237 PetscMPIInt size; 2238 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2239 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2240 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2241 PetscInt *masked,nmask,tmp,bs2,ishift; 2242 PetscScalar *aa; 2243 MPI_Comm comm = ((PetscObject)viewer)->comm; 2244 2245 PetscFunctionBegin; 2246 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2247 bs2 = bs*bs; 2248 2249 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2250 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2251 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2252 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2253 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2254 M = header[1]; N = header[2]; nz = header[3]; 2255 2256 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2257 2258 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2259 2260 /* 2261 This code adds extra rows to make sure the number of rows is 2262 divisible by the blocksize 2263 */ 2264 mbs = M/bs; 2265 extra_rows = bs - M + bs*(mbs); 2266 if (extra_rows == bs) extra_rows = 0; 2267 else mbs++; 2268 if (extra_rows) { 2269 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2270 } 2271 2272 /* Set global sizes if not already set */ 2273 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2274 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2275 } else { /* Check if the matrix global sizes are correct */ 2276 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2277 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); 2278 } 2279 2280 /* read in row lengths */ 2281 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2282 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2283 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2284 2285 /* read in column indices */ 2286 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2287 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2288 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2289 2290 /* loop over row lengths determining block row lengths */ 2291 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2292 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2293 ierr = PetscMalloc2(mbs,PetscInt,&mask,mbs,PetscInt,&masked);CHKERRQ(ierr); 2294 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2295 rowcount = 0; 2296 nzcount = 0; 2297 for (i=0; i<mbs; i++) { 2298 nmask = 0; 2299 for (j=0; j<bs; j++) { 2300 kmax = rowlengths[rowcount]; 2301 for (k=0; k<kmax; k++) { 2302 tmp = jj[nzcount++]/bs; /* block col. index */ 2303 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2304 } 2305 rowcount++; 2306 } 2307 s_browlengths[i] += nmask; 2308 2309 /* zero out the mask elements we set */ 2310 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2311 } 2312 2313 /* Do preallocation */ 2314 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2315 a = (Mat_SeqSBAIJ*)newmat->data; 2316 2317 /* set matrix "i" values */ 2318 a->i[0] = 0; 2319 for (i=1; i<= mbs; i++) { 2320 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2321 a->ilen[i-1] = s_browlengths[i-1]; 2322 } 2323 a->nz = a->i[mbs]; 2324 2325 /* read in nonzero values */ 2326 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2327 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2328 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2329 2330 /* set "a" and "j" values into matrix */ 2331 nzcount = 0; jcount = 0; 2332 for (i=0; i<mbs; i++) { 2333 nzcountb = nzcount; 2334 nmask = 0; 2335 for (j=0; j<bs; j++) { 2336 kmax = rowlengths[i*bs+j]; 2337 for (k=0; k<kmax; k++) { 2338 tmp = jj[nzcount++]/bs; /* block col. index */ 2339 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2340 } 2341 } 2342 /* sort the masked values */ 2343 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2344 2345 /* set "j" values into matrix */ 2346 maskcount = 1; 2347 for (j=0; j<nmask; j++) { 2348 a->j[jcount++] = masked[j]; 2349 mask[masked[j]] = maskcount++; 2350 } 2351 2352 /* set "a" values into matrix */ 2353 ishift = bs2*a->i[i]; 2354 for (j=0; j<bs; j++) { 2355 kmax = rowlengths[i*bs+j]; 2356 for (k=0; k<kmax; k++) { 2357 tmp = jj[nzcountb]/bs ; /* block col. index */ 2358 if (tmp >= i){ 2359 block = mask[tmp] - 1; 2360 point = jj[nzcountb] - bs*tmp; 2361 idx = ishift + bs2*block + j + bs*point; 2362 a->a[idx] = aa[nzcountb]; 2363 } 2364 nzcountb++; 2365 } 2366 } 2367 /* zero out the mask elements we set */ 2368 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2369 } 2370 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2371 2372 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2373 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2374 ierr = PetscFree(aa);CHKERRQ(ierr); 2375 ierr = PetscFree(jj);CHKERRQ(ierr); 2376 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2377 2378 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2379 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2380 ierr = MatView_Private(newmat);CHKERRQ(ierr); 2381 PetscFunctionReturn(0); 2382 } 2383 2384 #undef __FUNCT__ 2385 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2386 /*@ 2387 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2388 (upper triangular entries in CSR format) provided by the user. 2389 2390 Collective on MPI_Comm 2391 2392 Input Parameters: 2393 + comm - must be an MPI communicator of size 1 2394 . bs - size of block 2395 . m - number of rows 2396 . n - number of columns 2397 . i - row indices 2398 . j - column indices 2399 - a - matrix values 2400 2401 Output Parameter: 2402 . mat - the matrix 2403 2404 Level: advanced 2405 2406 Notes: 2407 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2408 once the matrix is destroyed 2409 2410 You cannot set new nonzero locations into this matrix, that will generate an error. 2411 2412 The i and j indices are 0 based 2413 2414 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 2415 it is the regular CSR format excluding the lower triangular elements. 2416 2417 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2418 2419 @*/ 2420 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2421 { 2422 PetscErrorCode ierr; 2423 PetscInt ii; 2424 Mat_SeqSBAIJ *sbaij; 2425 2426 PetscFunctionBegin; 2427 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2428 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2429 2430 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2431 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2432 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2433 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2434 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2435 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2436 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2437 2438 sbaij->i = i; 2439 sbaij->j = j; 2440 sbaij->a = a; 2441 sbaij->singlemalloc = PETSC_FALSE; 2442 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2443 sbaij->free_a = PETSC_FALSE; 2444 sbaij->free_ij = PETSC_FALSE; 2445 2446 for (ii=0; ii<m; ii++) { 2447 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2448 #if defined(PETSC_USE_DEBUG) 2449 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]); 2450 #endif 2451 } 2452 #if defined(PETSC_USE_DEBUG) 2453 for (ii=0; ii<sbaij->i[m]; ii++) { 2454 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2455 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]); 2456 } 2457 #endif 2458 2459 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2460 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2461 PetscFunctionReturn(0); 2462 } 2463 2464 2465 2466 2467 2468