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