1 /* 2 Defines the basic matrix operations for the AIJ (compressed row) 3 matrix storage format. 4 */ 5 6 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 #include <petscbt.h> 9 #include <petsc/private/kernels/blocktranspose.h> 10 11 PetscErrorCode MatSeqAIJSetTypeFromOptions(Mat A) 12 { 13 PetscBool flg; 14 char type[256]; 15 16 PetscFunctionBegin; 17 PetscObjectOptionsBegin((PetscObject)A); 18 PetscCall(PetscOptionsFList("-mat_seqaij_type","Matrix SeqAIJ type","MatSeqAIJSetType",MatSeqAIJList,"seqaij",type,256,&flg)); 19 if (flg) PetscCall(MatSeqAIJSetType(A,type)); 20 PetscOptionsEnd(); 21 PetscFunctionReturn(0); 22 } 23 24 PetscErrorCode MatGetColumnReductions_SeqAIJ(Mat A,PetscInt type,PetscReal *reductions) 25 { 26 PetscInt i,m,n; 27 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 28 29 PetscFunctionBegin; 30 PetscCall(MatGetSize(A,&m,&n)); 31 PetscCall(PetscArrayzero(reductions,n)); 32 if (type == NORM_2) { 33 for (i=0; i<aij->i[m]; i++) { 34 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 35 } 36 } else if (type == NORM_1) { 37 for (i=0; i<aij->i[m]; i++) { 38 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]); 39 } 40 } else if (type == NORM_INFINITY) { 41 for (i=0; i<aij->i[m]; i++) { 42 reductions[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),reductions[aij->j[i]]); 43 } 44 } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) { 45 for (i=0; i<aij->i[m]; i++) { 46 reductions[aij->j[i]] += PetscRealPart(aij->a[i]); 47 } 48 } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) { 49 for (i=0; i<aij->i[m]; i++) { 50 reductions[aij->j[i]] += PetscImaginaryPart(aij->a[i]); 51 } 52 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown reduction type"); 53 54 if (type == NORM_2) { 55 for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]); 56 } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) { 57 for (i=0; i<n; i++) reductions[i] /= m; 58 } 59 PetscFunctionReturn(0); 60 } 61 62 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is) 63 { 64 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 65 PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs; 66 const PetscInt *jj = a->j,*ii = a->i; 67 PetscInt *rows; 68 69 PetscFunctionBegin; 70 for (i=0; i<m; i++) { 71 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 72 cnt++; 73 } 74 } 75 PetscCall(PetscMalloc1(cnt,&rows)); 76 cnt = 0; 77 for (i=0; i<m; i++) { 78 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 79 rows[cnt] = i; 80 cnt++; 81 } 82 } 83 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is)); 84 PetscFunctionReturn(0); 85 } 86 87 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 88 { 89 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 90 const MatScalar *aa; 91 PetscInt i,m=A->rmap->n,cnt = 0; 92 const PetscInt *ii = a->i,*jj = a->j,*diag; 93 PetscInt *rows; 94 95 PetscFunctionBegin; 96 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 97 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 98 diag = a->diag; 99 for (i=0; i<m; i++) { 100 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 101 cnt++; 102 } 103 } 104 PetscCall(PetscMalloc1(cnt,&rows)); 105 cnt = 0; 106 for (i=0; i<m; i++) { 107 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 108 rows[cnt++] = i; 109 } 110 } 111 *nrows = cnt; 112 *zrows = rows; 113 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 118 { 119 PetscInt nrows,*rows; 120 121 PetscFunctionBegin; 122 *zrows = NULL; 123 PetscCall(MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows)); 124 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows)); 125 PetscFunctionReturn(0); 126 } 127 128 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 129 { 130 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 131 const MatScalar *aa; 132 PetscInt m=A->rmap->n,cnt = 0; 133 const PetscInt *ii; 134 PetscInt n,i,j,*rows; 135 136 PetscFunctionBegin; 137 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 138 *keptrows = NULL; 139 ii = a->i; 140 for (i=0; i<m; i++) { 141 n = ii[i+1] - ii[i]; 142 if (!n) { 143 cnt++; 144 goto ok1; 145 } 146 for (j=ii[i]; j<ii[i+1]; j++) { 147 if (aa[j] != 0.0) goto ok1; 148 } 149 cnt++; 150 ok1:; 151 } 152 if (!cnt) { 153 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 154 PetscFunctionReturn(0); 155 } 156 PetscCall(PetscMalloc1(A->rmap->n-cnt,&rows)); 157 cnt = 0; 158 for (i=0; i<m; i++) { 159 n = ii[i+1] - ii[i]; 160 if (!n) continue; 161 for (j=ii[i]; j<ii[i+1]; j++) { 162 if (aa[j] != 0.0) { 163 rows[cnt++] = i; 164 break; 165 } 166 } 167 } 168 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 169 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows)); 170 PetscFunctionReturn(0); 171 } 172 173 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 174 { 175 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 176 PetscInt i,m = Y->rmap->n; 177 const PetscInt *diag; 178 MatScalar *aa; 179 const PetscScalar *v; 180 PetscBool missing; 181 182 PetscFunctionBegin; 183 if (Y->assembled) { 184 PetscCall(MatMissingDiagonal_SeqAIJ(Y,&missing,NULL)); 185 if (!missing) { 186 diag = aij->diag; 187 PetscCall(VecGetArrayRead(D,&v)); 188 PetscCall(MatSeqAIJGetArray(Y,&aa)); 189 if (is == INSERT_VALUES) { 190 for (i=0; i<m; i++) { 191 aa[diag[i]] = v[i]; 192 } 193 } else { 194 for (i=0; i<m; i++) { 195 aa[diag[i]] += v[i]; 196 } 197 } 198 PetscCall(MatSeqAIJRestoreArray(Y,&aa)); 199 PetscCall(VecRestoreArrayRead(D,&v)); 200 PetscFunctionReturn(0); 201 } 202 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 203 } 204 PetscCall(MatDiagonalSet_Default(Y,D,is)); 205 PetscFunctionReturn(0); 206 } 207 208 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 209 { 210 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 211 PetscInt i,ishift; 212 213 PetscFunctionBegin; 214 *m = A->rmap->n; 215 if (!ia) PetscFunctionReturn(0); 216 ishift = 0; 217 if (symmetric && !A->structurally_symmetric) { 218 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja)); 219 } else if (oshift == 1) { 220 PetscInt *tia; 221 PetscInt nz = a->i[A->rmap->n]; 222 /* malloc space and add 1 to i and j indices */ 223 PetscCall(PetscMalloc1(A->rmap->n+1,&tia)); 224 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 225 *ia = tia; 226 if (ja) { 227 PetscInt *tja; 228 PetscCall(PetscMalloc1(nz+1,&tja)); 229 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 230 *ja = tja; 231 } 232 } else { 233 *ia = a->i; 234 if (ja) *ja = a->j; 235 } 236 PetscFunctionReturn(0); 237 } 238 239 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 240 { 241 PetscFunctionBegin; 242 if (!ia) PetscFunctionReturn(0); 243 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 244 PetscCall(PetscFree(*ia)); 245 if (ja) PetscCall(PetscFree(*ja)); 246 } 247 PetscFunctionReturn(0); 248 } 249 250 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 251 { 252 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 253 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 254 PetscInt nz = a->i[m],row,*jj,mr,col; 255 256 PetscFunctionBegin; 257 *nn = n; 258 if (!ia) PetscFunctionReturn(0); 259 if (symmetric) { 260 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,0,oshift,(PetscInt**)ia,(PetscInt**)ja)); 261 } else { 262 PetscCall(PetscCalloc1(n,&collengths)); 263 PetscCall(PetscMalloc1(n+1,&cia)); 264 PetscCall(PetscMalloc1(nz,&cja)); 265 jj = a->j; 266 for (i=0; i<nz; i++) { 267 collengths[jj[i]]++; 268 } 269 cia[0] = oshift; 270 for (i=0; i<n; i++) { 271 cia[i+1] = cia[i] + collengths[i]; 272 } 273 PetscCall(PetscArrayzero(collengths,n)); 274 jj = a->j; 275 for (row=0; row<m; row++) { 276 mr = a->i[row+1] - a->i[row]; 277 for (i=0; i<mr; i++) { 278 col = *jj++; 279 280 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 281 } 282 } 283 PetscCall(PetscFree(collengths)); 284 *ia = cia; *ja = cja; 285 } 286 PetscFunctionReturn(0); 287 } 288 289 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 290 { 291 PetscFunctionBegin; 292 if (!ia) PetscFunctionReturn(0); 293 294 PetscCall(PetscFree(*ia)); 295 PetscCall(PetscFree(*ja)); 296 PetscFunctionReturn(0); 297 } 298 299 /* 300 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 301 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 302 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 303 */ 304 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 305 { 306 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 307 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 308 PetscInt nz = a->i[m],row,mr,col,tmp; 309 PetscInt *cspidx; 310 const PetscInt *jj; 311 312 PetscFunctionBegin; 313 *nn = n; 314 if (!ia) PetscFunctionReturn(0); 315 316 PetscCall(PetscCalloc1(n,&collengths)); 317 PetscCall(PetscMalloc1(n+1,&cia)); 318 PetscCall(PetscMalloc1(nz,&cja)); 319 PetscCall(PetscMalloc1(nz,&cspidx)); 320 jj = a->j; 321 for (i=0; i<nz; i++) { 322 collengths[jj[i]]++; 323 } 324 cia[0] = oshift; 325 for (i=0; i<n; i++) { 326 cia[i+1] = cia[i] + collengths[i]; 327 } 328 PetscCall(PetscArrayzero(collengths,n)); 329 jj = a->j; 330 for (row=0; row<m; row++) { 331 mr = a->i[row+1] - a->i[row]; 332 for (i=0; i<mr; i++) { 333 col = *jj++; 334 tmp = cia[col] + collengths[col]++ - oshift; 335 cspidx[tmp] = a->i[row] + i; /* index of a->j */ 336 cja[tmp] = row + oshift; 337 } 338 } 339 PetscCall(PetscFree(collengths)); 340 *ia = cia; 341 *ja = cja; 342 *spidx = cspidx; 343 PetscFunctionReturn(0); 344 } 345 346 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 347 { 348 PetscFunctionBegin; 349 PetscCall(MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done)); 350 PetscCall(PetscFree(*spidx)); 351 PetscFunctionReturn(0); 352 } 353 354 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 355 { 356 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 357 PetscInt *ai = a->i; 358 PetscScalar *aa; 359 360 PetscFunctionBegin; 361 PetscCall(MatSeqAIJGetArray(A,&aa)); 362 PetscCall(PetscArraycpy(aa+ai[row],v,ai[row+1]-ai[row])); 363 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 364 PetscFunctionReturn(0); 365 } 366 367 /* 368 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 369 370 - a single row of values is set with each call 371 - no row or column indices are negative or (in error) larger than the number of rows or columns 372 - the values are always added to the matrix, not set 373 - no new locations are introduced in the nonzero structure of the matrix 374 375 This does NOT assume the global column indices are sorted 376 377 */ 378 379 #include <petsc/private/isimpl.h> 380 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 381 { 382 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 383 PetscInt low,high,t,row,nrow,i,col,l; 384 const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j; 385 PetscInt lastcol = -1; 386 MatScalar *ap,value,*aa; 387 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 388 389 PetscFunctionBegin; 390 PetscCall(MatSeqAIJGetArray(A,&aa)); 391 row = ridx[im[0]]; 392 rp = aj + ai[row]; 393 ap = aa + ai[row]; 394 nrow = ailen[row]; 395 low = 0; 396 high = nrow; 397 for (l=0; l<n; l++) { /* loop over added columns */ 398 col = cidx[in[l]]; 399 value = v[l]; 400 401 if (col <= lastcol) low = 0; 402 else high = nrow; 403 lastcol = col; 404 while (high-low > 5) { 405 t = (low+high)/2; 406 if (rp[t] > col) high = t; 407 else low = t; 408 } 409 for (i=low; i<high; i++) { 410 if (rp[i] == col) { 411 ap[i] += value; 412 low = i + 1; 413 break; 414 } 415 } 416 } 417 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 418 return 0; 419 } 420 421 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 422 { 423 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 424 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 425 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 426 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 427 MatScalar *ap=NULL,value=0.0,*aa; 428 PetscBool ignorezeroentries = a->ignorezeroentries; 429 PetscBool roworiented = a->roworiented; 430 431 PetscFunctionBegin; 432 PetscCall(MatSeqAIJGetArray(A,&aa)); 433 for (k=0; k<m; k++) { /* loop over added rows */ 434 row = im[k]; 435 if (row < 0) continue; 436 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 437 rp = aj + ai[row]; 438 if (!A->structure_only) ap = aa + ai[row]; 439 rmax = imax[row]; nrow = ailen[row]; 440 low = 0; 441 high = nrow; 442 for (l=0; l<n; l++) { /* loop over added columns */ 443 if (in[l] < 0) continue; 444 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 445 col = in[l]; 446 if (v && !A->structure_only) value = roworiented ? v[l + k*n] : v[k + l*m]; 447 if (!A->structure_only && value == 0.0 && ignorezeroentries && is == ADD_VALUES && row != col) continue; 448 449 if (col <= lastcol) low = 0; 450 else high = nrow; 451 lastcol = col; 452 while (high-low > 5) { 453 t = (low+high)/2; 454 if (rp[t] > col) high = t; 455 else low = t; 456 } 457 for (i=low; i<high; i++) { 458 if (rp[i] > col) break; 459 if (rp[i] == col) { 460 if (!A->structure_only) { 461 if (is == ADD_VALUES) { 462 ap[i] += value; 463 (void)PetscLogFlops(1.0); 464 } 465 else ap[i] = value; 466 } 467 low = i + 1; 468 goto noinsert; 469 } 470 } 471 if (value == 0.0 && ignorezeroentries && row != col) goto noinsert; 472 if (nonew == 1) goto noinsert; 473 PetscCheck(nonew != -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%" PetscInt_FMT ",%" PetscInt_FMT ") in the matrix",row,col); 474 if (A->structure_only) { 475 MatSeqXAIJReallocateAIJ_structure_only(A,A->rmap->n,1,nrow,row,col,rmax,ai,aj,rp,imax,nonew,MatScalar); 476 } else { 477 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 478 } 479 N = nrow++ - 1; a->nz++; high++; 480 /* shift up all the later entries in this row */ 481 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 482 rp[i] = col; 483 if (!A->structure_only) { 484 PetscCall(PetscArraymove(ap+i+1,ap+i,N-i+1)); 485 ap[i] = value; 486 } 487 low = i + 1; 488 A->nonzerostate++; 489 noinsert:; 490 } 491 ailen[row] = nrow; 492 } 493 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 494 PetscFunctionReturn(0); 495 } 496 497 PetscErrorCode MatSetValues_SeqAIJ_SortedFullNoPreallocation(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 498 { 499 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 500 PetscInt *rp,k,row; 501 PetscInt *ai = a->i; 502 PetscInt *aj = a->j; 503 MatScalar *aa,*ap; 504 505 PetscFunctionBegin; 506 PetscCheck(!A->was_assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot call on assembled matrix."); 507 PetscCheck(m*n+a->nz <= a->maxnz,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of entries in matrix will be larger than maximum nonzeros allocated for %" PetscInt_FMT " in MatSeqAIJSetTotalPreallocation()",a->maxnz); 508 509 PetscCall(MatSeqAIJGetArray(A,&aa)); 510 for (k=0; k<m; k++) { /* loop over added rows */ 511 row = im[k]; 512 rp = aj + ai[row]; 513 ap = aa + ai[row]; 514 515 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 516 if (!A->structure_only) { 517 if (v) { 518 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 519 v += n; 520 } else { 521 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 522 } 523 } 524 a->ilen[row] = n; 525 a->imax[row] = n; 526 a->i[row+1] = a->i[row]+n; 527 a->nz += n; 528 } 529 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 530 PetscFunctionReturn(0); 531 } 532 533 /*@ 534 MatSeqAIJSetTotalPreallocation - Sets an upper bound on the total number of expected nonzeros in the matrix. 535 536 Input Parameters: 537 + A - the SeqAIJ matrix 538 - nztotal - bound on the number of nonzeros 539 540 Level: advanced 541 542 Notes: 543 This can be called if you will be provided the matrix row by row (from row zero) with sorted column indices for each row. 544 Simply call MatSetValues() after this call to provide the matrix entries in the usual manner. This matrix may be used 545 as always with multiple matrix assemblies. 546 547 .seealso: `MatSetOption()`, `MAT_SORTED_FULL`, `MatSetValues()`, `MatSeqAIJSetPreallocation()` 548 @*/ 549 550 PetscErrorCode MatSeqAIJSetTotalPreallocation(Mat A,PetscInt nztotal) 551 { 552 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 553 554 PetscFunctionBegin; 555 PetscCall(PetscLayoutSetUp(A->rmap)); 556 PetscCall(PetscLayoutSetUp(A->cmap)); 557 a->maxnz = nztotal; 558 if (!a->imax) { 559 PetscCall(PetscMalloc1(A->rmap->n,&a->imax)); 560 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 561 } 562 if (!a->ilen) { 563 PetscCall(PetscMalloc1(A->rmap->n,&a->ilen)); 564 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 565 } else { 566 PetscCall(PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt))); 567 } 568 569 /* allocate the matrix space */ 570 if (A->structure_only) { 571 PetscCall(PetscMalloc1(nztotal,&a->j)); 572 PetscCall(PetscMalloc1(A->rmap->n+1,&a->i)); 573 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*sizeof(PetscInt))); 574 } else { 575 PetscCall(PetscMalloc3(nztotal,&a->a,nztotal,&a->j,A->rmap->n+1,&a->i)); 576 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*(sizeof(PetscScalar)+sizeof(PetscInt)))); 577 } 578 a->i[0] = 0; 579 if (A->structure_only) { 580 a->singlemalloc = PETSC_FALSE; 581 a->free_a = PETSC_FALSE; 582 } else { 583 a->singlemalloc = PETSC_TRUE; 584 a->free_a = PETSC_TRUE; 585 } 586 a->free_ij = PETSC_TRUE; 587 A->ops->setvalues = MatSetValues_SeqAIJ_SortedFullNoPreallocation; 588 A->preallocated = PETSC_TRUE; 589 PetscFunctionReturn(0); 590 } 591 592 PetscErrorCode MatSetValues_SeqAIJ_SortedFull(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 593 { 594 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 595 PetscInt *rp,k,row; 596 PetscInt *ai = a->i,*ailen = a->ilen; 597 PetscInt *aj = a->j; 598 MatScalar *aa,*ap; 599 600 PetscFunctionBegin; 601 PetscCall(MatSeqAIJGetArray(A,&aa)); 602 for (k=0; k<m; k++) { /* loop over added rows */ 603 row = im[k]; 604 PetscCheck(n <= a->imax[row],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Preallocation for row %" PetscInt_FMT " does not match number of columns provided",n); 605 rp = aj + ai[row]; 606 ap = aa + ai[row]; 607 if (!A->was_assembled) { 608 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 609 } 610 if (!A->structure_only) { 611 if (v) { 612 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 613 v += n; 614 } else { 615 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 616 } 617 } 618 ailen[row] = n; 619 a->nz += n; 620 } 621 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 622 PetscFunctionReturn(0); 623 } 624 625 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 626 { 627 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 628 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 629 PetscInt *ai = a->i,*ailen = a->ilen; 630 MatScalar *ap,*aa; 631 632 PetscFunctionBegin; 633 PetscCall(MatSeqAIJGetArray(A,&aa)); 634 for (k=0; k<m; k++) { /* loop over rows */ 635 row = im[k]; 636 if (row < 0) {v += n; continue;} /* negative row */ 637 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 638 rp = aj + ai[row]; ap = aa + ai[row]; 639 nrow = ailen[row]; 640 for (l=0; l<n; l++) { /* loop over columns */ 641 if (in[l] < 0) {v++; continue;} /* negative column */ 642 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 643 col = in[l]; 644 high = nrow; low = 0; /* assume unsorted */ 645 while (high-low > 5) { 646 t = (low+high)/2; 647 if (rp[t] > col) high = t; 648 else low = t; 649 } 650 for (i=low; i<high; i++) { 651 if (rp[i] > col) break; 652 if (rp[i] == col) { 653 *v++ = ap[i]; 654 goto finished; 655 } 656 } 657 *v++ = 0.0; 658 finished:; 659 } 660 } 661 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 662 PetscFunctionReturn(0); 663 } 664 665 PetscErrorCode MatView_SeqAIJ_Binary(Mat mat,PetscViewer viewer) 666 { 667 Mat_SeqAIJ *A = (Mat_SeqAIJ*)mat->data; 668 const PetscScalar *av; 669 PetscInt header[4],M,N,m,nz,i; 670 PetscInt *rowlens; 671 672 PetscFunctionBegin; 673 PetscCall(PetscViewerSetUp(viewer)); 674 675 M = mat->rmap->N; 676 N = mat->cmap->N; 677 m = mat->rmap->n; 678 nz = A->nz; 679 680 /* write matrix header */ 681 header[0] = MAT_FILE_CLASSID; 682 header[1] = M; header[2] = N; header[3] = nz; 683 PetscCall(PetscViewerBinaryWrite(viewer,header,4,PETSC_INT)); 684 685 /* fill in and store row lengths */ 686 PetscCall(PetscMalloc1(m,&rowlens)); 687 for (i=0; i<m; i++) rowlens[i] = A->i[i+1] - A->i[i]; 688 PetscCall(PetscViewerBinaryWrite(viewer,rowlens,m,PETSC_INT)); 689 PetscCall(PetscFree(rowlens)); 690 /* store column indices */ 691 PetscCall(PetscViewerBinaryWrite(viewer,A->j,nz,PETSC_INT)); 692 /* store nonzero values */ 693 PetscCall(MatSeqAIJGetArrayRead(mat,&av)); 694 PetscCall(PetscViewerBinaryWrite(viewer,av,nz,PETSC_SCALAR)); 695 PetscCall(MatSeqAIJRestoreArrayRead(mat,&av)); 696 697 /* write block size option to the viewer's .info file */ 698 PetscCall(MatView_Binary_BlockSizes(mat,viewer)); 699 PetscFunctionReturn(0); 700 } 701 702 static PetscErrorCode MatView_SeqAIJ_ASCII_structonly(Mat A,PetscViewer viewer) 703 { 704 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 705 PetscInt i,k,m=A->rmap->N; 706 707 PetscFunctionBegin; 708 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 709 for (i=0; i<m; i++) { 710 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 711 for (k=a->i[i]; k<a->i[i+1]; k++) { 712 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ") ",a->j[k])); 713 } 714 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 715 } 716 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 717 PetscFunctionReturn(0); 718 } 719 720 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 721 722 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 723 { 724 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 725 const PetscScalar *av; 726 PetscInt i,j,m = A->rmap->n; 727 const char *name; 728 PetscViewerFormat format; 729 730 PetscFunctionBegin; 731 if (A->structure_only) { 732 PetscCall(MatView_SeqAIJ_ASCII_structonly(A,viewer)); 733 PetscFunctionReturn(0); 734 } 735 736 PetscCall(PetscViewerGetFormat(viewer,&format)); 737 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 738 739 /* trigger copy to CPU if needed */ 740 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 741 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 742 if (format == PETSC_VIEWER_ASCII_MATLAB) { 743 PetscInt nofinalvalue = 0; 744 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) { 745 /* Need a dummy value to ensure the dimension of the matrix. */ 746 nofinalvalue = 1; 747 } 748 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 749 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n",m,A->cmap->n)); 750 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %" PetscInt_FMT " \n",a->nz)); 751 #if defined(PETSC_USE_COMPLEX) 752 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",4);\n",a->nz+nofinalvalue)); 753 #else 754 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",a->nz+nofinalvalue)); 755 #endif 756 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n")); 757 758 for (i=0; i<m; i++) { 759 for (j=a->i[i]; j<a->i[i+1]; j++) { 760 #if defined(PETSC_USE_COMPLEX) 761 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",i+1,a->j[j]+1,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 762 #else 763 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,a->j[j]+1,(double)a->a[j])); 764 #endif 765 } 766 } 767 if (nofinalvalue) { 768 #if defined(PETSC_USE_COMPLEX) 769 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",m,A->cmap->n,0.,0.)); 770 #else 771 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",m,A->cmap->n,0.0)); 772 #endif 773 } 774 PetscCall(PetscObjectGetName((PetscObject)A,&name)); 775 PetscCall(PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name)); 776 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 777 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 778 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 779 for (i=0; i<m; i++) { 780 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 781 for (j=a->i[i]; j<a->i[i+1]; j++) { 782 #if defined(PETSC_USE_COMPLEX) 783 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 784 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 785 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 786 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 787 } else if (PetscRealPart(a->a[j]) != 0.0) { 788 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 789 } 790 #else 791 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 792 #endif 793 } 794 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 795 } 796 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 797 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 798 PetscInt nzd=0,fshift=1,*sptr; 799 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 800 PetscCall(PetscMalloc1(m+1,&sptr)); 801 for (i=0; i<m; i++) { 802 sptr[i] = nzd+1; 803 for (j=a->i[i]; j<a->i[i+1]; j++) { 804 if (a->j[j] >= i) { 805 #if defined(PETSC_USE_COMPLEX) 806 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 807 #else 808 if (a->a[j] != 0.0) nzd++; 809 #endif 810 } 811 } 812 } 813 sptr[m] = nzd+1; 814 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n\n",m,nzd)); 815 for (i=0; i<m+1; i+=6) { 816 if (i+4<m) { 817 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5])); 818 } else if (i+3<m) { 819 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4])); 820 } else if (i+2<m) { 821 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3])); 822 } else if (i+1<m) { 823 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2])); 824 } else if (i<m) { 825 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1])); 826 } else { 827 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT "\n",sptr[i])); 828 } 829 } 830 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 831 PetscCall(PetscFree(sptr)); 832 for (i=0; i<m; i++) { 833 for (j=a->i[i]; j<a->i[i+1]; j++) { 834 if (a->j[j] >= i) PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " ",a->j[j]+fshift)); 835 } 836 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 837 } 838 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 839 for (i=0; i<m; i++) { 840 for (j=a->i[i]; j<a->i[i+1]; j++) { 841 if (a->j[j] >= i) { 842 #if defined(PETSC_USE_COMPLEX) 843 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 844 PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 845 } 846 #else 847 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j])); 848 #endif 849 } 850 } 851 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 852 } 853 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 854 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 855 PetscInt cnt = 0,jcnt; 856 PetscScalar value; 857 #if defined(PETSC_USE_COMPLEX) 858 PetscBool realonly = PETSC_TRUE; 859 860 for (i=0; i<a->i[m]; i++) { 861 if (PetscImaginaryPart(a->a[i]) != 0.0) { 862 realonly = PETSC_FALSE; 863 break; 864 } 865 } 866 #endif 867 868 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 869 for (i=0; i<m; i++) { 870 jcnt = 0; 871 for (j=0; j<A->cmap->n; j++) { 872 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 873 value = a->a[cnt++]; 874 jcnt++; 875 } else { 876 value = 0.0; 877 } 878 #if defined(PETSC_USE_COMPLEX) 879 if (realonly) { 880 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value))); 881 } else { 882 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value))); 883 } 884 #else 885 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value)); 886 #endif 887 } 888 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 889 } 890 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 891 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 892 PetscInt fshift=1; 893 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 894 #if defined(PETSC_USE_COMPLEX) 895 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n")); 896 #else 897 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n")); 898 #endif 899 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", m, A->cmap->n, a->nz)); 900 for (i=0; i<m; i++) { 901 for (j=a->i[i]; j<a->i[i+1]; j++) { 902 #if defined(PETSC_USE_COMPLEX) 903 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g %g\n", i+fshift,a->j[j]+fshift,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 904 #else 905 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j])); 906 #endif 907 } 908 } 909 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 910 } else { 911 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 912 if (A->factortype) { 913 for (i=0; i<m; i++) { 914 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 915 /* L part */ 916 for (j=a->i[i]; j<a->i[i+1]; j++) { 917 #if defined(PETSC_USE_COMPLEX) 918 if (PetscImaginaryPart(a->a[j]) > 0.0) { 919 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 920 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 921 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 922 } else { 923 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 924 } 925 #else 926 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 927 #endif 928 } 929 /* diagonal */ 930 j = a->diag[i]; 931 #if defined(PETSC_USE_COMPLEX) 932 if (PetscImaginaryPart(a->a[j]) > 0.0) { 933 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)PetscImaginaryPart(1.0/a->a[j]))); 934 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 935 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)(-PetscImaginaryPart(1.0/a->a[j])))); 936 } else { 937 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]))); 938 } 939 #else 940 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)(1.0/a->a[j]))); 941 #endif 942 943 /* U part */ 944 for (j=a->diag[i+1]+1; j<a->diag[i]; j++) { 945 #if defined(PETSC_USE_COMPLEX) 946 if (PetscImaginaryPart(a->a[j]) > 0.0) { 947 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 948 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 949 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 950 } else { 951 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 952 } 953 #else 954 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 955 #endif 956 } 957 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 958 } 959 } else { 960 for (i=0; i<m; i++) { 961 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 962 for (j=a->i[i]; j<a->i[i+1]; j++) { 963 #if defined(PETSC_USE_COMPLEX) 964 if (PetscImaginaryPart(a->a[j]) > 0.0) { 965 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 966 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 967 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 968 } else { 969 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 970 } 971 #else 972 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 973 #endif 974 } 975 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 976 } 977 } 978 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 979 } 980 PetscCall(PetscViewerFlush(viewer)); 981 PetscFunctionReturn(0); 982 } 983 984 #include <petscdraw.h> 985 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 986 { 987 Mat A = (Mat) Aa; 988 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 989 PetscInt i,j,m = A->rmap->n; 990 int color; 991 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 992 PetscViewer viewer; 993 PetscViewerFormat format; 994 const PetscScalar *aa; 995 996 PetscFunctionBegin; 997 PetscCall(PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer)); 998 PetscCall(PetscViewerGetFormat(viewer,&format)); 999 PetscCall(PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr)); 1000 1001 /* loop over matrix elements drawing boxes */ 1002 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1003 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1004 PetscDrawCollectiveBegin(draw); 1005 /* Blue for negative, Cyan for zero and Red for positive */ 1006 color = PETSC_DRAW_BLUE; 1007 for (i=0; i<m; i++) { 1008 y_l = m - i - 1.0; y_r = y_l + 1.0; 1009 for (j=a->i[i]; j<a->i[i+1]; j++) { 1010 x_l = a->j[j]; x_r = x_l + 1.0; 1011 if (PetscRealPart(aa[j]) >= 0.) continue; 1012 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1013 } 1014 } 1015 color = PETSC_DRAW_CYAN; 1016 for (i=0; i<m; i++) { 1017 y_l = m - i - 1.0; y_r = y_l + 1.0; 1018 for (j=a->i[i]; j<a->i[i+1]; j++) { 1019 x_l = a->j[j]; x_r = x_l + 1.0; 1020 if (aa[j] != 0.) continue; 1021 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1022 } 1023 } 1024 color = PETSC_DRAW_RED; 1025 for (i=0; i<m; i++) { 1026 y_l = m - i - 1.0; y_r = y_l + 1.0; 1027 for (j=a->i[i]; j<a->i[i+1]; j++) { 1028 x_l = a->j[j]; x_r = x_l + 1.0; 1029 if (PetscRealPart(aa[j]) <= 0.) continue; 1030 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1031 } 1032 } 1033 PetscDrawCollectiveEnd(draw); 1034 } else { 1035 /* use contour shading to indicate magnitude of values */ 1036 /* first determine max of all nonzero values */ 1037 PetscReal minv = 0.0, maxv = 0.0; 1038 PetscInt nz = a->nz, count = 0; 1039 PetscDraw popup; 1040 1041 for (i=0; i<nz; i++) { 1042 if (PetscAbsScalar(aa[i]) > maxv) maxv = PetscAbsScalar(aa[i]); 1043 } 1044 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1045 PetscCall(PetscDrawGetPopup(draw,&popup)); 1046 PetscCall(PetscDrawScalePopup(popup,minv,maxv)); 1047 1048 PetscDrawCollectiveBegin(draw); 1049 for (i=0; i<m; i++) { 1050 y_l = m - i - 1.0; 1051 y_r = y_l + 1.0; 1052 for (j=a->i[i]; j<a->i[i+1]; j++) { 1053 x_l = a->j[j]; 1054 x_r = x_l + 1.0; 1055 color = PetscDrawRealToColor(PetscAbsScalar(aa[count]),minv,maxv); 1056 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1057 count++; 1058 } 1059 } 1060 PetscDrawCollectiveEnd(draw); 1061 } 1062 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1063 PetscFunctionReturn(0); 1064 } 1065 1066 #include <petscdraw.h> 1067 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 1068 { 1069 PetscDraw draw; 1070 PetscReal xr,yr,xl,yl,h,w; 1071 PetscBool isnull; 1072 1073 PetscFunctionBegin; 1074 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 1075 PetscCall(PetscDrawIsNull(draw,&isnull)); 1076 if (isnull) PetscFunctionReturn(0); 1077 1078 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1079 xr += w; yr += h; xl = -w; yl = -h; 1080 PetscCall(PetscDrawSetCoordinates(draw,xl,yl,xr,yr)); 1081 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer)); 1082 PetscCall(PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A)); 1083 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL)); 1084 PetscCall(PetscDrawSave(draw)); 1085 PetscFunctionReturn(0); 1086 } 1087 1088 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 1089 { 1090 PetscBool iascii,isbinary,isdraw; 1091 1092 PetscFunctionBegin; 1093 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1094 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 1095 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 1096 if (iascii) { 1097 PetscCall(MatView_SeqAIJ_ASCII(A,viewer)); 1098 } else if (isbinary) { 1099 PetscCall(MatView_SeqAIJ_Binary(A,viewer)); 1100 } else if (isdraw) { 1101 PetscCall(MatView_SeqAIJ_Draw(A,viewer)); 1102 } 1103 PetscCall(MatView_SeqAIJ_Inode(A,viewer)); 1104 PetscFunctionReturn(0); 1105 } 1106 1107 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 1108 { 1109 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1110 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 1111 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 1112 MatScalar *aa = a->a,*ap; 1113 PetscReal ratio = 0.6; 1114 1115 PetscFunctionBegin; 1116 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 1117 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1118 if (A->was_assembled && A->ass_nonzerostate == A->nonzerostate) { 1119 /* we need to respect users asking to use or not the inodes routine in between matrix assemblies */ 1120 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1121 PetscFunctionReturn(0); 1122 } 1123 1124 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 1125 for (i=1; i<m; i++) { 1126 /* move each row back by the amount of empty slots (fshift) before it*/ 1127 fshift += imax[i-1] - ailen[i-1]; 1128 rmax = PetscMax(rmax,ailen[i]); 1129 if (fshift) { 1130 ip = aj + ai[i]; 1131 ap = aa + ai[i]; 1132 N = ailen[i]; 1133 PetscCall(PetscArraymove(ip-fshift,ip,N)); 1134 if (!A->structure_only) { 1135 PetscCall(PetscArraymove(ap-fshift,ap,N)); 1136 } 1137 } 1138 ai[i] = ai[i-1] + ailen[i-1]; 1139 } 1140 if (m) { 1141 fshift += imax[m-1] - ailen[m-1]; 1142 ai[m] = ai[m-1] + ailen[m-1]; 1143 } 1144 1145 /* reset ilen and imax for each row */ 1146 a->nonzerorowcnt = 0; 1147 if (A->structure_only) { 1148 PetscCall(PetscFree(a->imax)); 1149 PetscCall(PetscFree(a->ilen)); 1150 } else { /* !A->structure_only */ 1151 for (i=0; i<m; i++) { 1152 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1153 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 1154 } 1155 } 1156 a->nz = ai[m]; 1157 PetscCheck(!fshift || a->nounused != -1,PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, fshift); 1158 1159 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1160 PetscCall(PetscInfo(A,"Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded,%" PetscInt_FMT " used\n",m,A->cmap->n,fshift,a->nz)); 1161 PetscCall(PetscInfo(A,"Number of mallocs during MatSetValues() is %" PetscInt_FMT "\n",a->reallocs)); 1162 PetscCall(PetscInfo(A,"Maximum nonzeros in any row is %" PetscInt_FMT "\n",rmax)); 1163 1164 A->info.mallocs += a->reallocs; 1165 a->reallocs = 0; 1166 A->info.nz_unneeded = (PetscReal)fshift; 1167 a->rmax = rmax; 1168 1169 if (!A->structure_only) { 1170 PetscCall(MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio)); 1171 } 1172 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1173 PetscFunctionReturn(0); 1174 } 1175 1176 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1177 { 1178 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1179 PetscInt i,nz = a->nz; 1180 MatScalar *aa; 1181 1182 PetscFunctionBegin; 1183 PetscCall(MatSeqAIJGetArray(A,&aa)); 1184 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1185 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1186 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1187 PetscFunctionReturn(0); 1188 } 1189 1190 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1191 { 1192 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1193 PetscInt i,nz = a->nz; 1194 MatScalar *aa; 1195 1196 PetscFunctionBegin; 1197 PetscCall(MatSeqAIJGetArray(A,&aa)); 1198 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1199 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1200 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1201 PetscFunctionReturn(0); 1202 } 1203 1204 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1205 { 1206 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1207 MatScalar *aa; 1208 1209 PetscFunctionBegin; 1210 PetscCall(MatSeqAIJGetArrayWrite(A,&aa)); 1211 PetscCall(PetscArrayzero(aa,a->i[A->rmap->n])); 1212 PetscCall(MatSeqAIJRestoreArrayWrite(A,&aa)); 1213 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1214 PetscFunctionReturn(0); 1215 } 1216 1217 PETSC_INTERN PetscErrorCode MatResetPreallocationCOO_SeqAIJ(Mat A) 1218 { 1219 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1220 1221 PetscFunctionBegin; 1222 PetscCall(PetscFree(a->perm)); 1223 PetscCall(PetscFree(a->jmap)); 1224 PetscFunctionReturn(0); 1225 } 1226 1227 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1228 { 1229 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1230 1231 PetscFunctionBegin; 1232 #if defined(PETSC_USE_LOG) 1233 PetscLogObjectState((PetscObject)A,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT,A->rmap->n,A->cmap->n,a->nz); 1234 #endif 1235 PetscCall(MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i)); 1236 PetscCall(MatResetPreallocationCOO_SeqAIJ(A)); 1237 PetscCall(ISDestroy(&a->row)); 1238 PetscCall(ISDestroy(&a->col)); 1239 PetscCall(PetscFree(a->diag)); 1240 PetscCall(PetscFree(a->ibdiag)); 1241 PetscCall(PetscFree(a->imax)); 1242 PetscCall(PetscFree(a->ilen)); 1243 PetscCall(PetscFree(a->ipre)); 1244 PetscCall(PetscFree3(a->idiag,a->mdiag,a->ssor_work)); 1245 PetscCall(PetscFree(a->solve_work)); 1246 PetscCall(ISDestroy(&a->icol)); 1247 PetscCall(PetscFree(a->saved_values)); 1248 PetscCall(PetscFree2(a->compressedrow.i,a->compressedrow.rindex)); 1249 PetscCall(MatDestroy_SeqAIJ_Inode(A)); 1250 PetscCall(PetscFree(A->data)); 1251 1252 /* MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted may allocate this. 1253 That function is so heavily used (sometimes in an hidden way through multnumeric function pointers) 1254 that is hard to properly add this data to the MatProduct data. We free it here to avoid 1255 users reusing the matrix object with different data to incur in obscure segmentation faults 1256 due to different matrix sizes */ 1257 PetscCall(PetscObjectCompose((PetscObject)A,"__PETSc__ab_dense",NULL)); 1258 1259 PetscCall(PetscObjectChangeTypeName((PetscObject)A,NULL)); 1260 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL)); 1261 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL)); 1262 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL)); 1263 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL)); 1264 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL)); 1265 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL)); 1266 #if defined(PETSC_HAVE_CUDA) 1267 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcusparse_C",NULL)); 1268 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",NULL)); 1269 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",NULL)); 1270 #endif 1271 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 1272 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijkokkos_C",NULL)); 1273 #endif 1274 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcrl_C",NULL)); 1275 #if defined(PETSC_HAVE_ELEMENTAL) 1276 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL)); 1277 #endif 1278 #if defined(PETSC_HAVE_SCALAPACK) 1279 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_scalapack_C",NULL)); 1280 #endif 1281 #if defined(PETSC_HAVE_HYPRE) 1282 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_hypre_C",NULL)); 1283 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",NULL)); 1284 #endif 1285 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL)); 1286 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsell_C",NULL)); 1287 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_is_C",NULL)); 1288 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL)); 1289 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL)); 1290 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatResetPreallocation_C",NULL)); 1291 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL)); 1292 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL)); 1293 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_is_seqaij_C",NULL)); 1294 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqdense_seqaij_C",NULL)); 1295 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaij_C",NULL)); 1296 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJKron_C",NULL)); 1297 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL)); 1298 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL)); 1299 PetscFunctionReturn(0); 1300 } 1301 1302 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1303 { 1304 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1305 1306 PetscFunctionBegin; 1307 switch (op) { 1308 case MAT_ROW_ORIENTED: 1309 a->roworiented = flg; 1310 break; 1311 case MAT_KEEP_NONZERO_PATTERN: 1312 a->keepnonzeropattern = flg; 1313 break; 1314 case MAT_NEW_NONZERO_LOCATIONS: 1315 a->nonew = (flg ? 0 : 1); 1316 break; 1317 case MAT_NEW_NONZERO_LOCATION_ERR: 1318 a->nonew = (flg ? -1 : 0); 1319 break; 1320 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1321 a->nonew = (flg ? -2 : 0); 1322 break; 1323 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1324 a->nounused = (flg ? -1 : 0); 1325 break; 1326 case MAT_IGNORE_ZERO_ENTRIES: 1327 a->ignorezeroentries = flg; 1328 break; 1329 case MAT_SPD: 1330 case MAT_SYMMETRIC: 1331 case MAT_STRUCTURALLY_SYMMETRIC: 1332 case MAT_HERMITIAN: 1333 case MAT_SYMMETRY_ETERNAL: 1334 case MAT_STRUCTURE_ONLY: 1335 /* These options are handled directly by MatSetOption() */ 1336 break; 1337 case MAT_FORCE_DIAGONAL_ENTRIES: 1338 case MAT_IGNORE_OFF_PROC_ENTRIES: 1339 case MAT_USE_HASH_TABLE: 1340 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 1341 break; 1342 case MAT_USE_INODES: 1343 PetscCall(MatSetOption_SeqAIJ_Inode(A,MAT_USE_INODES,flg)); 1344 break; 1345 case MAT_SUBMAT_SINGLEIS: 1346 A->submat_singleis = flg; 1347 break; 1348 case MAT_SORTED_FULL: 1349 if (flg) A->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 1350 else A->ops->setvalues = MatSetValues_SeqAIJ; 1351 break; 1352 case MAT_FORM_EXPLICIT_TRANSPOSE: 1353 A->form_explicit_transpose = flg; 1354 break; 1355 default: 1356 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1357 } 1358 PetscFunctionReturn(0); 1359 } 1360 1361 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1362 { 1363 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1364 PetscInt i,j,n,*ai=a->i,*aj=a->j; 1365 PetscScalar *x; 1366 const PetscScalar *aa; 1367 1368 PetscFunctionBegin; 1369 PetscCall(VecGetLocalSize(v,&n)); 1370 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1371 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1372 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1373 PetscInt *diag=a->diag; 1374 PetscCall(VecGetArrayWrite(v,&x)); 1375 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1376 PetscCall(VecRestoreArrayWrite(v,&x)); 1377 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1378 PetscFunctionReturn(0); 1379 } 1380 1381 PetscCall(VecGetArrayWrite(v,&x)); 1382 for (i=0; i<n; i++) { 1383 x[i] = 0.0; 1384 for (j=ai[i]; j<ai[i+1]; j++) { 1385 if (aj[j] == i) { 1386 x[i] = aa[j]; 1387 break; 1388 } 1389 } 1390 } 1391 PetscCall(VecRestoreArrayWrite(v,&x)); 1392 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1393 PetscFunctionReturn(0); 1394 } 1395 1396 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1397 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1398 { 1399 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1400 const MatScalar *aa; 1401 PetscScalar *y; 1402 const PetscScalar *x; 1403 PetscInt m = A->rmap->n; 1404 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1405 const MatScalar *v; 1406 PetscScalar alpha; 1407 PetscInt n,i,j; 1408 const PetscInt *idx,*ii,*ridx=NULL; 1409 Mat_CompressedRow cprow = a->compressedrow; 1410 PetscBool usecprow = cprow.use; 1411 #endif 1412 1413 PetscFunctionBegin; 1414 if (zz != yy) PetscCall(VecCopy(zz,yy)); 1415 PetscCall(VecGetArrayRead(xx,&x)); 1416 PetscCall(VecGetArray(yy,&y)); 1417 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1418 1419 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1420 fortranmulttransposeaddaij_(&m,x,a->i,a->j,aa,y); 1421 #else 1422 if (usecprow) { 1423 m = cprow.nrows; 1424 ii = cprow.i; 1425 ridx = cprow.rindex; 1426 } else { 1427 ii = a->i; 1428 } 1429 for (i=0; i<m; i++) { 1430 idx = a->j + ii[i]; 1431 v = aa + ii[i]; 1432 n = ii[i+1] - ii[i]; 1433 if (usecprow) { 1434 alpha = x[ridx[i]]; 1435 } else { 1436 alpha = x[i]; 1437 } 1438 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1439 } 1440 #endif 1441 PetscCall(PetscLogFlops(2.0*a->nz)); 1442 PetscCall(VecRestoreArrayRead(xx,&x)); 1443 PetscCall(VecRestoreArray(yy,&y)); 1444 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1445 PetscFunctionReturn(0); 1446 } 1447 1448 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1449 { 1450 PetscFunctionBegin; 1451 PetscCall(VecSet(yy,0.0)); 1452 PetscCall(MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy)); 1453 PetscFunctionReturn(0); 1454 } 1455 1456 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1457 1458 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1459 { 1460 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1461 PetscScalar *y; 1462 const PetscScalar *x; 1463 const MatScalar *aa,*a_a; 1464 PetscInt m=A->rmap->n; 1465 const PetscInt *aj,*ii,*ridx=NULL; 1466 PetscInt n,i; 1467 PetscScalar sum; 1468 PetscBool usecprow=a->compressedrow.use; 1469 1470 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1471 #pragma disjoint(*x,*y,*aa) 1472 #endif 1473 1474 PetscFunctionBegin; 1475 if (a->inode.use && a->inode.checked) { 1476 PetscCall(MatMult_SeqAIJ_Inode(A,xx,yy)); 1477 PetscFunctionReturn(0); 1478 } 1479 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1480 PetscCall(VecGetArrayRead(xx,&x)); 1481 PetscCall(VecGetArray(yy,&y)); 1482 ii = a->i; 1483 if (usecprow) { /* use compressed row format */ 1484 PetscCall(PetscArrayzero(y,m)); 1485 m = a->compressedrow.nrows; 1486 ii = a->compressedrow.i; 1487 ridx = a->compressedrow.rindex; 1488 for (i=0; i<m; i++) { 1489 n = ii[i+1] - ii[i]; 1490 aj = a->j + ii[i]; 1491 aa = a_a + ii[i]; 1492 sum = 0.0; 1493 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1494 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1495 y[*ridx++] = sum; 1496 } 1497 } else { /* do not use compressed row format */ 1498 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1499 aj = a->j; 1500 aa = a_a; 1501 fortranmultaij_(&m,x,ii,aj,aa,y); 1502 #else 1503 for (i=0; i<m; i++) { 1504 n = ii[i+1] - ii[i]; 1505 aj = a->j + ii[i]; 1506 aa = a_a + ii[i]; 1507 sum = 0.0; 1508 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1509 y[i] = sum; 1510 } 1511 #endif 1512 } 1513 PetscCall(PetscLogFlops(2.0*a->nz - a->nonzerorowcnt)); 1514 PetscCall(VecRestoreArrayRead(xx,&x)); 1515 PetscCall(VecRestoreArray(yy,&y)); 1516 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1517 PetscFunctionReturn(0); 1518 } 1519 1520 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1521 { 1522 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1523 PetscScalar *y; 1524 const PetscScalar *x; 1525 const MatScalar *aa,*a_a; 1526 PetscInt m=A->rmap->n; 1527 const PetscInt *aj,*ii,*ridx=NULL; 1528 PetscInt n,i,nonzerorow=0; 1529 PetscScalar sum; 1530 PetscBool usecprow=a->compressedrow.use; 1531 1532 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1533 #pragma disjoint(*x,*y,*aa) 1534 #endif 1535 1536 PetscFunctionBegin; 1537 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1538 PetscCall(VecGetArrayRead(xx,&x)); 1539 PetscCall(VecGetArray(yy,&y)); 1540 if (usecprow) { /* use compressed row format */ 1541 m = a->compressedrow.nrows; 1542 ii = a->compressedrow.i; 1543 ridx = a->compressedrow.rindex; 1544 for (i=0; i<m; i++) { 1545 n = ii[i+1] - ii[i]; 1546 aj = a->j + ii[i]; 1547 aa = a_a + ii[i]; 1548 sum = 0.0; 1549 nonzerorow += (n>0); 1550 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1551 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1552 y[*ridx++] = sum; 1553 } 1554 } else { /* do not use compressed row format */ 1555 ii = a->i; 1556 for (i=0; i<m; i++) { 1557 n = ii[i+1] - ii[i]; 1558 aj = a->j + ii[i]; 1559 aa = a_a + ii[i]; 1560 sum = 0.0; 1561 nonzerorow += (n>0); 1562 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1563 y[i] = sum; 1564 } 1565 } 1566 PetscCall(PetscLogFlops(2.0*a->nz - nonzerorow)); 1567 PetscCall(VecRestoreArrayRead(xx,&x)); 1568 PetscCall(VecRestoreArray(yy,&y)); 1569 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1570 PetscFunctionReturn(0); 1571 } 1572 1573 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1574 { 1575 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1576 PetscScalar *y,*z; 1577 const PetscScalar *x; 1578 const MatScalar *aa,*a_a; 1579 PetscInt m = A->rmap->n,*aj,*ii; 1580 PetscInt n,i,*ridx=NULL; 1581 PetscScalar sum; 1582 PetscBool usecprow=a->compressedrow.use; 1583 1584 PetscFunctionBegin; 1585 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1586 PetscCall(VecGetArrayRead(xx,&x)); 1587 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1588 if (usecprow) { /* use compressed row format */ 1589 if (zz != yy) { 1590 PetscCall(PetscArraycpy(z,y,m)); 1591 } 1592 m = a->compressedrow.nrows; 1593 ii = a->compressedrow.i; 1594 ridx = a->compressedrow.rindex; 1595 for (i=0; i<m; i++) { 1596 n = ii[i+1] - ii[i]; 1597 aj = a->j + ii[i]; 1598 aa = a_a + ii[i]; 1599 sum = y[*ridx]; 1600 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1601 z[*ridx++] = sum; 1602 } 1603 } else { /* do not use compressed row format */ 1604 ii = a->i; 1605 for (i=0; i<m; i++) { 1606 n = ii[i+1] - ii[i]; 1607 aj = a->j + ii[i]; 1608 aa = a_a + ii[i]; 1609 sum = y[i]; 1610 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1611 z[i] = sum; 1612 } 1613 } 1614 PetscCall(PetscLogFlops(2.0*a->nz)); 1615 PetscCall(VecRestoreArrayRead(xx,&x)); 1616 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1617 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1618 PetscFunctionReturn(0); 1619 } 1620 1621 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1622 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1623 { 1624 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1625 PetscScalar *y,*z; 1626 const PetscScalar *x; 1627 const MatScalar *aa,*a_a; 1628 const PetscInt *aj,*ii,*ridx=NULL; 1629 PetscInt m = A->rmap->n,n,i; 1630 PetscScalar sum; 1631 PetscBool usecprow=a->compressedrow.use; 1632 1633 PetscFunctionBegin; 1634 if (a->inode.use && a->inode.checked) { 1635 PetscCall(MatMultAdd_SeqAIJ_Inode(A,xx,yy,zz)); 1636 PetscFunctionReturn(0); 1637 } 1638 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1639 PetscCall(VecGetArrayRead(xx,&x)); 1640 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1641 if (usecprow) { /* use compressed row format */ 1642 if (zz != yy) { 1643 PetscCall(PetscArraycpy(z,y,m)); 1644 } 1645 m = a->compressedrow.nrows; 1646 ii = a->compressedrow.i; 1647 ridx = a->compressedrow.rindex; 1648 for (i=0; i<m; i++) { 1649 n = ii[i+1] - ii[i]; 1650 aj = a->j + ii[i]; 1651 aa = a_a + ii[i]; 1652 sum = y[*ridx]; 1653 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1654 z[*ridx++] = sum; 1655 } 1656 } else { /* do not use compressed row format */ 1657 ii = a->i; 1658 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1659 aj = a->j; 1660 aa = a_a; 1661 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1662 #else 1663 for (i=0; i<m; i++) { 1664 n = ii[i+1] - ii[i]; 1665 aj = a->j + ii[i]; 1666 aa = a_a + ii[i]; 1667 sum = y[i]; 1668 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1669 z[i] = sum; 1670 } 1671 #endif 1672 } 1673 PetscCall(PetscLogFlops(2.0*a->nz)); 1674 PetscCall(VecRestoreArrayRead(xx,&x)); 1675 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1676 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1677 PetscFunctionReturn(0); 1678 } 1679 1680 /* 1681 Adds diagonal pointers to sparse matrix structure. 1682 */ 1683 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1684 { 1685 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1686 PetscInt i,j,m = A->rmap->n; 1687 PetscBool alreadySet = PETSC_TRUE; 1688 1689 PetscFunctionBegin; 1690 if (!a->diag) { 1691 PetscCall(PetscMalloc1(m,&a->diag)); 1692 PetscCall(PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt))); 1693 alreadySet = PETSC_FALSE; 1694 } 1695 for (i=0; i<A->rmap->n; i++) { 1696 /* If A's diagonal is already correctly set, this fast track enables cheap and repeated MatMarkDiagonal_SeqAIJ() calls */ 1697 if (alreadySet) { 1698 PetscInt pos = a->diag[i]; 1699 if (pos >= a->i[i] && pos < a->i[i+1] && a->j[pos] == i) continue; 1700 } 1701 1702 a->diag[i] = a->i[i+1]; 1703 for (j=a->i[i]; j<a->i[i+1]; j++) { 1704 if (a->j[j] == i) { 1705 a->diag[i] = j; 1706 break; 1707 } 1708 } 1709 } 1710 PetscFunctionReturn(0); 1711 } 1712 1713 PetscErrorCode MatShift_SeqAIJ(Mat A,PetscScalar v) 1714 { 1715 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1716 const PetscInt *diag = (const PetscInt*)a->diag; 1717 const PetscInt *ii = (const PetscInt*) a->i; 1718 PetscInt i,*mdiag = NULL; 1719 PetscInt cnt = 0; /* how many diagonals are missing */ 1720 1721 PetscFunctionBegin; 1722 if (!A->preallocated || !a->nz) { 1723 PetscCall(MatSeqAIJSetPreallocation(A,1,NULL)); 1724 PetscCall(MatShift_Basic(A,v)); 1725 PetscFunctionReturn(0); 1726 } 1727 1728 if (a->diagonaldense) { 1729 cnt = 0; 1730 } else { 1731 PetscCall(PetscCalloc1(A->rmap->n,&mdiag)); 1732 for (i=0; i<A->rmap->n; i++) { 1733 if (i < A->cmap->n && diag[i] >= ii[i+1]) { /* 'out of range' rows never have diagonals */ 1734 cnt++; 1735 mdiag[i] = 1; 1736 } 1737 } 1738 } 1739 if (!cnt) { 1740 PetscCall(MatShift_Basic(A,v)); 1741 } else { 1742 PetscScalar *olda = a->a; /* preserve pointers to current matrix nonzeros structure and values */ 1743 PetscInt *oldj = a->j, *oldi = a->i; 1744 PetscBool singlemalloc = a->singlemalloc,free_a = a->free_a,free_ij = a->free_ij; 1745 1746 a->a = NULL; 1747 a->j = NULL; 1748 a->i = NULL; 1749 /* increase the values in imax for each row where a diagonal is being inserted then reallocate the matrix data structures */ 1750 for (i=0; i<PetscMin(A->rmap->n,A->cmap->n); i++) { 1751 a->imax[i] += mdiag[i]; 1752 } 1753 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax)); 1754 1755 /* copy old values into new matrix data structure */ 1756 for (i=0; i<A->rmap->n; i++) { 1757 PetscCall(MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES)); 1758 if (i < A->cmap->n) { 1759 PetscCall(MatSetValue(A,i,i,v,ADD_VALUES)); 1760 } 1761 } 1762 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1763 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1764 if (singlemalloc) { 1765 PetscCall(PetscFree3(olda,oldj,oldi)); 1766 } else { 1767 if (free_a) PetscCall(PetscFree(olda)); 1768 if (free_ij) PetscCall(PetscFree(oldj)); 1769 if (free_ij) PetscCall(PetscFree(oldi)); 1770 } 1771 } 1772 PetscCall(PetscFree(mdiag)); 1773 a->diagonaldense = PETSC_TRUE; 1774 PetscFunctionReturn(0); 1775 } 1776 1777 /* 1778 Checks for missing diagonals 1779 */ 1780 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1781 { 1782 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1783 PetscInt *diag,*ii = a->i,i; 1784 1785 PetscFunctionBegin; 1786 *missing = PETSC_FALSE; 1787 if (A->rmap->n > 0 && !ii) { 1788 *missing = PETSC_TRUE; 1789 if (d) *d = 0; 1790 PetscCall(PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n")); 1791 } else { 1792 PetscInt n; 1793 n = PetscMin(A->rmap->n, A->cmap->n); 1794 diag = a->diag; 1795 for (i=0; i<n; i++) { 1796 if (diag[i] >= ii[i+1]) { 1797 *missing = PETSC_TRUE; 1798 if (d) *d = i; 1799 PetscCall(PetscInfo(A,"Matrix is missing diagonal number %" PetscInt_FMT "\n",i)); 1800 break; 1801 } 1802 } 1803 } 1804 PetscFunctionReturn(0); 1805 } 1806 1807 #include <petscblaslapack.h> 1808 #include <petsc/private/kernels/blockinvert.h> 1809 1810 /* 1811 Note that values is allocated externally by the PC and then passed into this routine 1812 */ 1813 PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag) 1814 { 1815 PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots; 1816 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 1817 const PetscReal shift = 0.0; 1818 PetscInt ipvt[5]; 1819 PetscScalar work[25],*v_work; 1820 1821 PetscFunctionBegin; 1822 allowzeropivot = PetscNot(A->erroriffailure); 1823 for (i=0; i<nblocks; i++) ncnt += bsizes[i]; 1824 PetscCheck(ncnt == n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Total blocksizes %" PetscInt_FMT " doesn't match number matrix rows %" PetscInt_FMT,ncnt,n); 1825 for (i=0; i<nblocks; i++) { 1826 bsizemax = PetscMax(bsizemax,bsizes[i]); 1827 } 1828 PetscCall(PetscMalloc1(bsizemax,&indx)); 1829 if (bsizemax > 7) { 1830 PetscCall(PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots)); 1831 } 1832 ncnt = 0; 1833 for (i=0; i<nblocks; i++) { 1834 for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j; 1835 PetscCall(MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag)); 1836 switch (bsizes[i]) { 1837 case 1: 1838 *diag = 1.0/(*diag); 1839 break; 1840 case 2: 1841 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 1842 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1843 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 1844 break; 1845 case 3: 1846 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 1847 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1848 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 1849 break; 1850 case 4: 1851 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 1852 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1853 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 1854 break; 1855 case 5: 1856 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 1857 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1858 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 1859 break; 1860 case 6: 1861 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 1862 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1863 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 1864 break; 1865 case 7: 1866 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 1867 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1868 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 1869 break; 1870 default: 1871 PetscCall(PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 1872 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1873 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bsizes[i])); 1874 } 1875 ncnt += bsizes[i]; 1876 diag += bsizes[i]*bsizes[i]; 1877 } 1878 if (bsizemax > 7) { 1879 PetscCall(PetscFree2(v_work,v_pivots)); 1880 } 1881 PetscCall(PetscFree(indx)); 1882 PetscFunctionReturn(0); 1883 } 1884 1885 /* 1886 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1887 */ 1888 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1889 { 1890 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1891 PetscInt i,*diag,m = A->rmap->n; 1892 const MatScalar *v; 1893 PetscScalar *idiag,*mdiag; 1894 1895 PetscFunctionBegin; 1896 if (a->idiagvalid) PetscFunctionReturn(0); 1897 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1898 diag = a->diag; 1899 if (!a->idiag) { 1900 PetscCall(PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work)); 1901 PetscCall(PetscLogObjectMemory((PetscObject)A,3*m*sizeof(PetscScalar))); 1902 } 1903 1904 mdiag = a->mdiag; 1905 idiag = a->idiag; 1906 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 1907 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1908 for (i=0; i<m; i++) { 1909 mdiag[i] = v[diag[i]]; 1910 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1911 if (PetscRealPart(fshift)) { 1912 PetscCall(PetscInfo(A,"Zero diagonal on row %" PetscInt_FMT "\n",i)); 1913 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1914 A->factorerror_zeropivot_value = 0.0; 1915 A->factorerror_zeropivot_row = i; 1916 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %" PetscInt_FMT,i); 1917 } 1918 idiag[i] = 1.0/v[diag[i]]; 1919 } 1920 PetscCall(PetscLogFlops(m)); 1921 } else { 1922 for (i=0; i<m; i++) { 1923 mdiag[i] = v[diag[i]]; 1924 idiag[i] = omega/(fshift + v[diag[i]]); 1925 } 1926 PetscCall(PetscLogFlops(2.0*m)); 1927 } 1928 a->idiagvalid = PETSC_TRUE; 1929 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 1930 PetscFunctionReturn(0); 1931 } 1932 1933 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1934 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1935 { 1936 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1937 PetscScalar *x,d,sum,*t,scale; 1938 const MatScalar *v,*idiag=NULL,*mdiag,*aa; 1939 const PetscScalar *b, *bs,*xb, *ts; 1940 PetscInt n,m = A->rmap->n,i; 1941 const PetscInt *idx,*diag; 1942 1943 PetscFunctionBegin; 1944 if (a->inode.use && a->inode.checked && omega == 1.0 && fshift == 0.0) { 1945 PetscCall(MatSOR_SeqAIJ_Inode(A,bb,omega,flag,fshift,its,lits,xx)); 1946 PetscFunctionReturn(0); 1947 } 1948 its = its*lits; 1949 1950 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1951 if (!a->idiagvalid) PetscCall(MatInvertDiagonal_SeqAIJ(A,omega,fshift)); 1952 a->fshift = fshift; 1953 a->omega = omega; 1954 1955 diag = a->diag; 1956 t = a->ssor_work; 1957 idiag = a->idiag; 1958 mdiag = a->mdiag; 1959 1960 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1961 PetscCall(VecGetArray(xx,&x)); 1962 PetscCall(VecGetArrayRead(bb,&b)); 1963 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1964 if (flag == SOR_APPLY_UPPER) { 1965 /* apply (U + D/omega) to the vector */ 1966 bs = b; 1967 for (i=0; i<m; i++) { 1968 d = fshift + mdiag[i]; 1969 n = a->i[i+1] - diag[i] - 1; 1970 idx = a->j + diag[i] + 1; 1971 v = aa + diag[i] + 1; 1972 sum = b[i]*d/omega; 1973 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1974 x[i] = sum; 1975 } 1976 PetscCall(VecRestoreArray(xx,&x)); 1977 PetscCall(VecRestoreArrayRead(bb,&b)); 1978 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1979 PetscCall(PetscLogFlops(a->nz)); 1980 PetscFunctionReturn(0); 1981 } 1982 1983 PetscCheck(flag != SOR_APPLY_LOWER,PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1984 else if (flag & SOR_EISENSTAT) { 1985 /* Let A = L + U + D; where L is lower triangular, 1986 U is upper triangular, E = D/omega; This routine applies 1987 1988 (L + E)^{-1} A (U + E)^{-1} 1989 1990 to a vector efficiently using Eisenstat's trick. 1991 */ 1992 scale = (2.0/omega) - 1.0; 1993 1994 /* x = (E + U)^{-1} b */ 1995 for (i=m-1; i>=0; i--) { 1996 n = a->i[i+1] - diag[i] - 1; 1997 idx = a->j + diag[i] + 1; 1998 v = aa + diag[i] + 1; 1999 sum = b[i]; 2000 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2001 x[i] = sum*idiag[i]; 2002 } 2003 2004 /* t = b - (2*E - D)x */ 2005 v = aa; 2006 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 2007 2008 /* t = (E + L)^{-1}t */ 2009 ts = t; 2010 diag = a->diag; 2011 for (i=0; i<m; i++) { 2012 n = diag[i] - a->i[i]; 2013 idx = a->j + a->i[i]; 2014 v = aa + a->i[i]; 2015 sum = t[i]; 2016 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 2017 t[i] = sum*idiag[i]; 2018 /* x = x + t */ 2019 x[i] += t[i]; 2020 } 2021 2022 PetscCall(PetscLogFlops(6.0*m-1 + 2.0*a->nz)); 2023 PetscCall(VecRestoreArray(xx,&x)); 2024 PetscCall(VecRestoreArrayRead(bb,&b)); 2025 PetscFunctionReturn(0); 2026 } 2027 if (flag & SOR_ZERO_INITIAL_GUESS) { 2028 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2029 for (i=0; i<m; i++) { 2030 n = diag[i] - a->i[i]; 2031 idx = a->j + a->i[i]; 2032 v = aa + a->i[i]; 2033 sum = b[i]; 2034 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2035 t[i] = sum; 2036 x[i] = sum*idiag[i]; 2037 } 2038 xb = t; 2039 PetscCall(PetscLogFlops(a->nz)); 2040 } else xb = b; 2041 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2042 for (i=m-1; i>=0; i--) { 2043 n = a->i[i+1] - diag[i] - 1; 2044 idx = a->j + diag[i] + 1; 2045 v = aa + diag[i] + 1; 2046 sum = xb[i]; 2047 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2048 if (xb == b) { 2049 x[i] = sum*idiag[i]; 2050 } else { 2051 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2052 } 2053 } 2054 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2055 } 2056 its--; 2057 } 2058 while (its--) { 2059 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2060 for (i=0; i<m; i++) { 2061 /* lower */ 2062 n = diag[i] - a->i[i]; 2063 idx = a->j + a->i[i]; 2064 v = aa + a->i[i]; 2065 sum = b[i]; 2066 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2067 t[i] = sum; /* save application of the lower-triangular part */ 2068 /* upper */ 2069 n = a->i[i+1] - diag[i] - 1; 2070 idx = a->j + diag[i] + 1; 2071 v = aa + diag[i] + 1; 2072 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2073 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2074 } 2075 xb = t; 2076 PetscCall(PetscLogFlops(2.0*a->nz)); 2077 } else xb = b; 2078 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2079 for (i=m-1; i>=0; i--) { 2080 sum = xb[i]; 2081 if (xb == b) { 2082 /* whole matrix (no checkpointing available) */ 2083 n = a->i[i+1] - a->i[i]; 2084 idx = a->j + a->i[i]; 2085 v = aa + a->i[i]; 2086 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2087 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 2088 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 2089 n = a->i[i+1] - diag[i] - 1; 2090 idx = a->j + diag[i] + 1; 2091 v = aa + diag[i] + 1; 2092 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2093 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2094 } 2095 } 2096 if (xb == b) { 2097 PetscCall(PetscLogFlops(2.0*a->nz)); 2098 } else { 2099 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2100 } 2101 } 2102 } 2103 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2104 PetscCall(VecRestoreArray(xx,&x)); 2105 PetscCall(VecRestoreArrayRead(bb,&b)); 2106 PetscFunctionReturn(0); 2107 } 2108 2109 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 2110 { 2111 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2112 2113 PetscFunctionBegin; 2114 info->block_size = 1.0; 2115 info->nz_allocated = a->maxnz; 2116 info->nz_used = a->nz; 2117 info->nz_unneeded = (a->maxnz - a->nz); 2118 info->assemblies = A->num_ass; 2119 info->mallocs = A->info.mallocs; 2120 info->memory = ((PetscObject)A)->mem; 2121 if (A->factortype) { 2122 info->fill_ratio_given = A->info.fill_ratio_given; 2123 info->fill_ratio_needed = A->info.fill_ratio_needed; 2124 info->factor_mallocs = A->info.factor_mallocs; 2125 } else { 2126 info->fill_ratio_given = 0; 2127 info->fill_ratio_needed = 0; 2128 info->factor_mallocs = 0; 2129 } 2130 PetscFunctionReturn(0); 2131 } 2132 2133 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2134 { 2135 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2136 PetscInt i,m = A->rmap->n - 1; 2137 const PetscScalar *xx; 2138 PetscScalar *bb,*aa; 2139 PetscInt d = 0; 2140 2141 PetscFunctionBegin; 2142 if (x && b) { 2143 PetscCall(VecGetArrayRead(x,&xx)); 2144 PetscCall(VecGetArray(b,&bb)); 2145 for (i=0; i<N; i++) { 2146 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2147 if (rows[i] >= A->cmap->n) continue; 2148 bb[rows[i]] = diag*xx[rows[i]]; 2149 } 2150 PetscCall(VecRestoreArrayRead(x,&xx)); 2151 PetscCall(VecRestoreArray(b,&bb)); 2152 } 2153 2154 PetscCall(MatSeqAIJGetArray(A,&aa)); 2155 if (a->keepnonzeropattern) { 2156 for (i=0; i<N; i++) { 2157 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2158 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2159 } 2160 if (diag != 0.0) { 2161 for (i=0; i<N; i++) { 2162 d = rows[i]; 2163 if (rows[i] >= A->cmap->n) continue; 2164 PetscCheck(a->diag[d] < a->i[d+1],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in the zeroed row %" PetscInt_FMT,d); 2165 } 2166 for (i=0; i<N; i++) { 2167 if (rows[i] >= A->cmap->n) continue; 2168 aa[a->diag[rows[i]]] = diag; 2169 } 2170 } 2171 } else { 2172 if (diag != 0.0) { 2173 for (i=0; i<N; i++) { 2174 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2175 if (a->ilen[rows[i]] > 0) { 2176 if (rows[i] >= A->cmap->n) { 2177 a->ilen[rows[i]] = 0; 2178 } else { 2179 a->ilen[rows[i]] = 1; 2180 aa[a->i[rows[i]]] = diag; 2181 a->j[a->i[rows[i]]] = rows[i]; 2182 } 2183 } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */ 2184 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2185 } 2186 } 2187 } else { 2188 for (i=0; i<N; i++) { 2189 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2190 a->ilen[rows[i]] = 0; 2191 } 2192 } 2193 A->nonzerostate++; 2194 } 2195 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2196 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2197 PetscFunctionReturn(0); 2198 } 2199 2200 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2201 { 2202 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2203 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2204 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2205 const PetscScalar *xx; 2206 PetscScalar *bb,*aa; 2207 2208 PetscFunctionBegin; 2209 if (!N) PetscFunctionReturn(0); 2210 PetscCall(MatSeqAIJGetArray(A,&aa)); 2211 if (x && b) { 2212 PetscCall(VecGetArrayRead(x,&xx)); 2213 PetscCall(VecGetArray(b,&bb)); 2214 vecs = PETSC_TRUE; 2215 } 2216 PetscCall(PetscCalloc1(A->rmap->n,&zeroed)); 2217 for (i=0; i<N; i++) { 2218 PetscCheck(rows[i] >= 0 && rows[i] <= m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2219 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2220 2221 zeroed[rows[i]] = PETSC_TRUE; 2222 } 2223 for (i=0; i<A->rmap->n; i++) { 2224 if (!zeroed[i]) { 2225 for (j=a->i[i]; j<a->i[i+1]; j++) { 2226 if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) { 2227 if (vecs) bb[i] -= aa[j]*xx[a->j[j]]; 2228 aa[j] = 0.0; 2229 } 2230 } 2231 } else if (vecs && i < A->cmap->N) bb[i] = diag*xx[i]; 2232 } 2233 if (x && b) { 2234 PetscCall(VecRestoreArrayRead(x,&xx)); 2235 PetscCall(VecRestoreArray(b,&bb)); 2236 } 2237 PetscCall(PetscFree(zeroed)); 2238 if (diag != 0.0) { 2239 PetscCall(MatMissingDiagonal_SeqAIJ(A,&missing,&d)); 2240 if (missing) { 2241 for (i=0; i<N; i++) { 2242 if (rows[i] >= A->cmap->N) continue; 2243 PetscCheck(!a->nonew || rows[i] < d,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %" PetscInt_FMT " (%" PetscInt_FMT ")",d,rows[i]); 2244 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2245 } 2246 } else { 2247 for (i=0; i<N; i++) { 2248 aa[a->diag[rows[i]]] = diag; 2249 } 2250 } 2251 } 2252 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2253 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2254 PetscFunctionReturn(0); 2255 } 2256 2257 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2258 { 2259 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2260 const PetscScalar *aa; 2261 PetscInt *itmp; 2262 2263 PetscFunctionBegin; 2264 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2265 *nz = a->i[row+1] - a->i[row]; 2266 if (v) *v = (PetscScalar*)(aa + a->i[row]); 2267 if (idx) { 2268 itmp = a->j + a->i[row]; 2269 if (*nz) *idx = itmp; 2270 else *idx = NULL; 2271 } 2272 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2273 PetscFunctionReturn(0); 2274 } 2275 2276 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2277 { 2278 PetscFunctionBegin; 2279 if (nz) *nz = 0; 2280 if (idx) *idx = NULL; 2281 if (v) *v = NULL; 2282 PetscFunctionReturn(0); 2283 } 2284 2285 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2286 { 2287 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2288 const MatScalar *v; 2289 PetscReal sum = 0.0; 2290 PetscInt i,j; 2291 2292 PetscFunctionBegin; 2293 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 2294 if (type == NORM_FROBENIUS) { 2295 #if defined(PETSC_USE_REAL___FP16) 2296 PetscBLASInt one = 1,nz = a->nz; 2297 PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&nz,v,&one)); 2298 #else 2299 for (i=0; i<a->nz; i++) { 2300 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2301 } 2302 *nrm = PetscSqrtReal(sum); 2303 #endif 2304 PetscCall(PetscLogFlops(2.0*a->nz)); 2305 } else if (type == NORM_1) { 2306 PetscReal *tmp; 2307 PetscInt *jj = a->j; 2308 PetscCall(PetscCalloc1(A->cmap->n+1,&tmp)); 2309 *nrm = 0.0; 2310 for (j=0; j<a->nz; j++) { 2311 tmp[*jj++] += PetscAbsScalar(*v); v++; 2312 } 2313 for (j=0; j<A->cmap->n; j++) { 2314 if (tmp[j] > *nrm) *nrm = tmp[j]; 2315 } 2316 PetscCall(PetscFree(tmp)); 2317 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2318 } else if (type == NORM_INFINITY) { 2319 *nrm = 0.0; 2320 for (j=0; j<A->rmap->n; j++) { 2321 const PetscScalar *v2 = v + a->i[j]; 2322 sum = 0.0; 2323 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2324 sum += PetscAbsScalar(*v2); v2++; 2325 } 2326 if (sum > *nrm) *nrm = sum; 2327 } 2328 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2329 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2330 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 2331 PetscFunctionReturn(0); 2332 } 2333 2334 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2335 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2336 { 2337 PetscInt i,j,anzj; 2338 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2339 PetscInt an=A->cmap->N,am=A->rmap->N; 2340 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2341 2342 PetscFunctionBegin; 2343 /* Allocate space for symbolic transpose info and work array */ 2344 PetscCall(PetscCalloc1(an+1,&ati)); 2345 PetscCall(PetscMalloc1(ai[am],&atj)); 2346 PetscCall(PetscMalloc1(an,&atfill)); 2347 2348 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2349 /* Note: offset by 1 for fast conversion into csr format. */ 2350 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2351 /* Form ati for csr format of A^T. */ 2352 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2353 2354 /* Copy ati into atfill so we have locations of the next free space in atj */ 2355 PetscCall(PetscArraycpy(atfill,ati,an)); 2356 2357 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2358 for (i=0;i<am;i++) { 2359 anzj = ai[i+1] - ai[i]; 2360 for (j=0;j<anzj;j++) { 2361 atj[atfill[*aj]] = i; 2362 atfill[*aj++] += 1; 2363 } 2364 } 2365 2366 /* Clean up temporary space and complete requests. */ 2367 PetscCall(PetscFree(atfill)); 2368 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B)); 2369 PetscCall(MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs))); 2370 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 2371 2372 b = (Mat_SeqAIJ*)((*B)->data); 2373 b->free_a = PETSC_FALSE; 2374 b->free_ij = PETSC_TRUE; 2375 b->nonew = 0; 2376 PetscFunctionReturn(0); 2377 } 2378 2379 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2380 { 2381 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2382 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2383 const MatScalar *va,*vb; 2384 PetscInt ma,na,mb,nb, i; 2385 2386 PetscFunctionBegin; 2387 PetscCall(MatGetSize(A,&ma,&na)); 2388 PetscCall(MatGetSize(B,&mb,&nb)); 2389 if (ma!=nb || na!=mb) { 2390 *f = PETSC_FALSE; 2391 PetscFunctionReturn(0); 2392 } 2393 PetscCall(MatSeqAIJGetArrayRead(A,&va)); 2394 PetscCall(MatSeqAIJGetArrayRead(B,&vb)); 2395 aii = aij->i; bii = bij->i; 2396 adx = aij->j; bdx = bij->j; 2397 PetscCall(PetscMalloc1(ma,&aptr)); 2398 PetscCall(PetscMalloc1(mb,&bptr)); 2399 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2400 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2401 2402 *f = PETSC_TRUE; 2403 for (i=0; i<ma; i++) { 2404 while (aptr[i]<aii[i+1]) { 2405 PetscInt idc,idr; 2406 PetscScalar vc,vr; 2407 /* column/row index/value */ 2408 idc = adx[aptr[i]]; 2409 idr = bdx[bptr[idc]]; 2410 vc = va[aptr[i]]; 2411 vr = vb[bptr[idc]]; 2412 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2413 *f = PETSC_FALSE; 2414 goto done; 2415 } else { 2416 aptr[i]++; 2417 if (B || i!=idc) bptr[idc]++; 2418 } 2419 } 2420 } 2421 done: 2422 PetscCall(PetscFree(aptr)); 2423 PetscCall(PetscFree(bptr)); 2424 PetscCall(MatSeqAIJRestoreArrayRead(A,&va)); 2425 PetscCall(MatSeqAIJRestoreArrayRead(B,&vb)); 2426 PetscFunctionReturn(0); 2427 } 2428 2429 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2430 { 2431 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2432 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2433 MatScalar *va,*vb; 2434 PetscInt ma,na,mb,nb, i; 2435 2436 PetscFunctionBegin; 2437 PetscCall(MatGetSize(A,&ma,&na)); 2438 PetscCall(MatGetSize(B,&mb,&nb)); 2439 if (ma!=nb || na!=mb) { 2440 *f = PETSC_FALSE; 2441 PetscFunctionReturn(0); 2442 } 2443 aii = aij->i; bii = bij->i; 2444 adx = aij->j; bdx = bij->j; 2445 va = aij->a; vb = bij->a; 2446 PetscCall(PetscMalloc1(ma,&aptr)); 2447 PetscCall(PetscMalloc1(mb,&bptr)); 2448 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2449 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2450 2451 *f = PETSC_TRUE; 2452 for (i=0; i<ma; i++) { 2453 while (aptr[i]<aii[i+1]) { 2454 PetscInt idc,idr; 2455 PetscScalar vc,vr; 2456 /* column/row index/value */ 2457 idc = adx[aptr[i]]; 2458 idr = bdx[bptr[idc]]; 2459 vc = va[aptr[i]]; 2460 vr = vb[bptr[idc]]; 2461 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2462 *f = PETSC_FALSE; 2463 goto done; 2464 } else { 2465 aptr[i]++; 2466 if (B || i!=idc) bptr[idc]++; 2467 } 2468 } 2469 } 2470 done: 2471 PetscCall(PetscFree(aptr)); 2472 PetscCall(PetscFree(bptr)); 2473 PetscFunctionReturn(0); 2474 } 2475 2476 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2477 { 2478 PetscFunctionBegin; 2479 PetscCall(MatIsTranspose_SeqAIJ(A,A,tol,f)); 2480 PetscFunctionReturn(0); 2481 } 2482 2483 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2484 { 2485 PetscFunctionBegin; 2486 PetscCall(MatIsHermitianTranspose_SeqAIJ(A,A,tol,f)); 2487 PetscFunctionReturn(0); 2488 } 2489 2490 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2491 { 2492 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2493 const PetscScalar *l,*r; 2494 PetscScalar x; 2495 MatScalar *v; 2496 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2497 const PetscInt *jj; 2498 2499 PetscFunctionBegin; 2500 if (ll) { 2501 /* The local size is used so that VecMPI can be passed to this routine 2502 by MatDiagonalScale_MPIAIJ */ 2503 PetscCall(VecGetLocalSize(ll,&m)); 2504 PetscCheck(m == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2505 PetscCall(VecGetArrayRead(ll,&l)); 2506 PetscCall(MatSeqAIJGetArray(A,&v)); 2507 for (i=0; i<m; i++) { 2508 x = l[i]; 2509 M = a->i[i+1] - a->i[i]; 2510 for (j=0; j<M; j++) (*v++) *= x; 2511 } 2512 PetscCall(VecRestoreArrayRead(ll,&l)); 2513 PetscCall(PetscLogFlops(nz)); 2514 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2515 } 2516 if (rr) { 2517 PetscCall(VecGetLocalSize(rr,&n)); 2518 PetscCheck(n == A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2519 PetscCall(VecGetArrayRead(rr,&r)); 2520 PetscCall(MatSeqAIJGetArray(A,&v)); 2521 jj = a->j; 2522 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2523 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2524 PetscCall(VecRestoreArrayRead(rr,&r)); 2525 PetscCall(PetscLogFlops(nz)); 2526 } 2527 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 2528 PetscFunctionReturn(0); 2529 } 2530 2531 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2532 { 2533 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2534 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2535 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2536 const PetscInt *irow,*icol; 2537 const PetscScalar *aa; 2538 PetscInt nrows,ncols; 2539 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2540 MatScalar *a_new,*mat_a; 2541 Mat C; 2542 PetscBool stride; 2543 2544 PetscFunctionBegin; 2545 PetscCall(ISGetIndices(isrow,&irow)); 2546 PetscCall(ISGetLocalSize(isrow,&nrows)); 2547 PetscCall(ISGetLocalSize(iscol,&ncols)); 2548 2549 PetscCall(PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride)); 2550 if (stride) { 2551 PetscCall(ISStrideGetInfo(iscol,&first,&step)); 2552 } else { 2553 first = 0; 2554 step = 0; 2555 } 2556 if (stride && step == 1) { 2557 /* special case of contiguous rows */ 2558 PetscCall(PetscMalloc2(nrows,&lens,nrows,&starts)); 2559 /* loop over new rows determining lens and starting points */ 2560 for (i=0; i<nrows; i++) { 2561 kstart = ai[irow[i]]; 2562 kend = kstart + ailen[irow[i]]; 2563 starts[i] = kstart; 2564 for (k=kstart; k<kend; k++) { 2565 if (aj[k] >= first) { 2566 starts[i] = k; 2567 break; 2568 } 2569 } 2570 sum = 0; 2571 while (k < kend) { 2572 if (aj[k++] >= first+ncols) break; 2573 sum++; 2574 } 2575 lens[i] = sum; 2576 } 2577 /* create submatrix */ 2578 if (scall == MAT_REUSE_MATRIX) { 2579 PetscInt n_cols,n_rows; 2580 PetscCall(MatGetSize(*B,&n_rows,&n_cols)); 2581 PetscCheck(n_rows == nrows && n_cols == ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2582 PetscCall(MatZeroEntries(*B)); 2583 C = *B; 2584 } else { 2585 PetscInt rbs,cbs; 2586 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2587 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2588 PetscCall(ISGetBlockSize(isrow,&rbs)); 2589 PetscCall(ISGetBlockSize(iscol,&cbs)); 2590 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2591 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2592 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2593 } 2594 c = (Mat_SeqAIJ*)C->data; 2595 2596 /* loop over rows inserting into submatrix */ 2597 a_new = c->a; 2598 j_new = c->j; 2599 i_new = c->i; 2600 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2601 for (i=0; i<nrows; i++) { 2602 ii = starts[i]; 2603 lensi = lens[i]; 2604 for (k=0; k<lensi; k++) { 2605 *j_new++ = aj[ii+k] - first; 2606 } 2607 PetscCall(PetscArraycpy(a_new,aa + starts[i],lensi)); 2608 a_new += lensi; 2609 i_new[i+1] = i_new[i] + lensi; 2610 c->ilen[i] = lensi; 2611 } 2612 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2613 PetscCall(PetscFree2(lens,starts)); 2614 } else { 2615 PetscCall(ISGetIndices(iscol,&icol)); 2616 PetscCall(PetscCalloc1(oldcols,&smap)); 2617 PetscCall(PetscMalloc1(1+nrows,&lens)); 2618 for (i=0; i<ncols; i++) { 2619 PetscCheck(icol[i] < oldcols,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requesting column beyond largest column icol[%" PetscInt_FMT "] %" PetscInt_FMT " >= A->cmap->n %" PetscInt_FMT,i,icol[i],oldcols); 2620 smap[icol[i]] = i+1; 2621 } 2622 2623 /* determine lens of each row */ 2624 for (i=0; i<nrows; i++) { 2625 kstart = ai[irow[i]]; 2626 kend = kstart + a->ilen[irow[i]]; 2627 lens[i] = 0; 2628 for (k=kstart; k<kend; k++) { 2629 if (smap[aj[k]]) { 2630 lens[i]++; 2631 } 2632 } 2633 } 2634 /* Create and fill new matrix */ 2635 if (scall == MAT_REUSE_MATRIX) { 2636 PetscBool equal; 2637 2638 c = (Mat_SeqAIJ*)((*B)->data); 2639 PetscCheck((*B)->rmap->n == nrows && (*B)->cmap->n == ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2640 PetscCall(PetscArraycmp(c->ilen,lens,(*B)->rmap->n,&equal)); 2641 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2642 PetscCall(PetscArrayzero(c->ilen,(*B)->rmap->n)); 2643 C = *B; 2644 } else { 2645 PetscInt rbs,cbs; 2646 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2647 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2648 PetscCall(ISGetBlockSize(isrow,&rbs)); 2649 PetscCall(ISGetBlockSize(iscol,&cbs)); 2650 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2651 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2652 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2653 } 2654 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2655 c = (Mat_SeqAIJ*)(C->data); 2656 for (i=0; i<nrows; i++) { 2657 row = irow[i]; 2658 kstart = ai[row]; 2659 kend = kstart + a->ilen[row]; 2660 mat_i = c->i[i]; 2661 mat_j = c->j + mat_i; 2662 mat_a = c->a + mat_i; 2663 mat_ilen = c->ilen + i; 2664 for (k=kstart; k<kend; k++) { 2665 if ((tcol=smap[a->j[k]])) { 2666 *mat_j++ = tcol - 1; 2667 *mat_a++ = aa[k]; 2668 (*mat_ilen)++; 2669 2670 } 2671 } 2672 } 2673 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2674 /* Free work space */ 2675 PetscCall(ISRestoreIndices(iscol,&icol)); 2676 PetscCall(PetscFree(smap)); 2677 PetscCall(PetscFree(lens)); 2678 /* sort */ 2679 for (i = 0; i < nrows; i++) { 2680 PetscInt ilen; 2681 2682 mat_i = c->i[i]; 2683 mat_j = c->j + mat_i; 2684 mat_a = c->a + mat_i; 2685 ilen = c->ilen[i]; 2686 PetscCall(PetscSortIntWithScalarArray(ilen,mat_j,mat_a)); 2687 } 2688 } 2689 #if defined(PETSC_HAVE_DEVICE) 2690 PetscCall(MatBindToCPU(C,A->boundtocpu)); 2691 #endif 2692 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 2693 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 2694 2695 PetscCall(ISRestoreIndices(isrow,&irow)); 2696 *B = C; 2697 PetscFunctionReturn(0); 2698 } 2699 2700 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2701 { 2702 Mat B; 2703 2704 PetscFunctionBegin; 2705 if (scall == MAT_INITIAL_MATRIX) { 2706 PetscCall(MatCreate(subComm,&B)); 2707 PetscCall(MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n)); 2708 PetscCall(MatSetBlockSizesFromMats(B,mat,mat)); 2709 PetscCall(MatSetType(B,MATSEQAIJ)); 2710 PetscCall(MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE)); 2711 *subMat = B; 2712 } else { 2713 PetscCall(MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN)); 2714 } 2715 PetscFunctionReturn(0); 2716 } 2717 2718 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2719 { 2720 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2721 Mat outA; 2722 PetscBool row_identity,col_identity; 2723 2724 PetscFunctionBegin; 2725 PetscCheck(info->levels == 0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2726 2727 PetscCall(ISIdentity(row,&row_identity)); 2728 PetscCall(ISIdentity(col,&col_identity)); 2729 2730 outA = inA; 2731 outA->factortype = MAT_FACTOR_LU; 2732 PetscCall(PetscFree(inA->solvertype)); 2733 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype)); 2734 2735 PetscCall(PetscObjectReference((PetscObject)row)); 2736 PetscCall(ISDestroy(&a->row)); 2737 2738 a->row = row; 2739 2740 PetscCall(PetscObjectReference((PetscObject)col)); 2741 PetscCall(ISDestroy(&a->col)); 2742 2743 a->col = col; 2744 2745 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2746 PetscCall(ISDestroy(&a->icol)); 2747 PetscCall(ISInvertPermutation(col,PETSC_DECIDE,&a->icol)); 2748 PetscCall(PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol)); 2749 2750 if (!a->solve_work) { /* this matrix may have been factored before */ 2751 PetscCall(PetscMalloc1(inA->rmap->n+1,&a->solve_work)); 2752 PetscCall(PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar))); 2753 } 2754 2755 PetscCall(MatMarkDiagonal_SeqAIJ(inA)); 2756 if (row_identity && col_identity) { 2757 PetscCall(MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info)); 2758 } else { 2759 PetscCall(MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info)); 2760 } 2761 PetscFunctionReturn(0); 2762 } 2763 2764 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2765 { 2766 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2767 PetscScalar *v; 2768 PetscBLASInt one = 1,bnz; 2769 2770 PetscFunctionBegin; 2771 PetscCall(MatSeqAIJGetArray(inA,&v)); 2772 PetscCall(PetscBLASIntCast(a->nz,&bnz)); 2773 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&alpha,v,&one)); 2774 PetscCall(PetscLogFlops(a->nz)); 2775 PetscCall(MatSeqAIJRestoreArray(inA,&v)); 2776 PetscCall(MatSeqAIJInvalidateDiagonal(inA)); 2777 PetscFunctionReturn(0); 2778 } 2779 2780 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2781 { 2782 PetscInt i; 2783 2784 PetscFunctionBegin; 2785 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2786 PetscCall(PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr)); 2787 2788 for (i=0; i<submatj->nrqr; ++i) { 2789 PetscCall(PetscFree(submatj->sbuf2[i])); 2790 } 2791 PetscCall(PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1)); 2792 2793 if (submatj->rbuf1) { 2794 PetscCall(PetscFree(submatj->rbuf1[0])); 2795 PetscCall(PetscFree(submatj->rbuf1)); 2796 } 2797 2798 for (i=0; i<submatj->nrqs; ++i) { 2799 PetscCall(PetscFree(submatj->rbuf3[i])); 2800 } 2801 PetscCall(PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3)); 2802 PetscCall(PetscFree(submatj->pa)); 2803 } 2804 2805 #if defined(PETSC_USE_CTABLE) 2806 PetscCall(PetscTableDestroy((PetscTable*)&submatj->rmap)); 2807 if (submatj->cmap_loc) PetscCall(PetscFree(submatj->cmap_loc)); 2808 PetscCall(PetscFree(submatj->rmap_loc)); 2809 #else 2810 PetscCall(PetscFree(submatj->rmap)); 2811 #endif 2812 2813 if (!submatj->allcolumns) { 2814 #if defined(PETSC_USE_CTABLE) 2815 PetscCall(PetscTableDestroy((PetscTable*)&submatj->cmap)); 2816 #else 2817 PetscCall(PetscFree(submatj->cmap)); 2818 #endif 2819 } 2820 PetscCall(PetscFree(submatj->row2proc)); 2821 2822 PetscCall(PetscFree(submatj)); 2823 PetscFunctionReturn(0); 2824 } 2825 2826 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2827 { 2828 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2829 Mat_SubSppt *submatj = c->submatis1; 2830 2831 PetscFunctionBegin; 2832 PetscCall((*submatj->destroy)(C)); 2833 PetscCall(MatDestroySubMatrix_Private(submatj)); 2834 PetscFunctionReturn(0); 2835 } 2836 2837 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2838 { 2839 PetscInt i; 2840 Mat C; 2841 Mat_SeqAIJ *c; 2842 Mat_SubSppt *submatj; 2843 2844 PetscFunctionBegin; 2845 for (i=0; i<n; i++) { 2846 C = (*mat)[i]; 2847 c = (Mat_SeqAIJ*)C->data; 2848 submatj = c->submatis1; 2849 if (submatj) { 2850 if (--((PetscObject)C)->refct <= 0) { 2851 PetscCall((*submatj->destroy)(C)); 2852 PetscCall(MatDestroySubMatrix_Private(submatj)); 2853 PetscCall(PetscFree(C->defaultvectype)); 2854 PetscCall(PetscLayoutDestroy(&C->rmap)); 2855 PetscCall(PetscLayoutDestroy(&C->cmap)); 2856 PetscCall(PetscHeaderDestroy(&C)); 2857 } 2858 } else { 2859 PetscCall(MatDestroy(&C)); 2860 } 2861 } 2862 2863 /* Destroy Dummy submatrices created for reuse */ 2864 PetscCall(MatDestroySubMatrices_Dummy(n,mat)); 2865 2866 PetscCall(PetscFree(*mat)); 2867 PetscFunctionReturn(0); 2868 } 2869 2870 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2871 { 2872 PetscInt i; 2873 2874 PetscFunctionBegin; 2875 if (scall == MAT_INITIAL_MATRIX) { 2876 PetscCall(PetscCalloc1(n+1,B)); 2877 } 2878 2879 for (i=0; i<n; i++) { 2880 PetscCall(MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i])); 2881 } 2882 PetscFunctionReturn(0); 2883 } 2884 2885 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2886 { 2887 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2888 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2889 const PetscInt *idx; 2890 PetscInt start,end,*ai,*aj; 2891 PetscBT table; 2892 2893 PetscFunctionBegin; 2894 m = A->rmap->n; 2895 ai = a->i; 2896 aj = a->j; 2897 2898 PetscCheck(ov >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2899 2900 PetscCall(PetscMalloc1(m+1,&nidx)); 2901 PetscCall(PetscBTCreate(m,&table)); 2902 2903 for (i=0; i<is_max; i++) { 2904 /* Initialize the two local arrays */ 2905 isz = 0; 2906 PetscCall(PetscBTMemzero(m,table)); 2907 2908 /* Extract the indices, assume there can be duplicate entries */ 2909 PetscCall(ISGetIndices(is[i],&idx)); 2910 PetscCall(ISGetLocalSize(is[i],&n)); 2911 2912 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2913 for (j=0; j<n; ++j) { 2914 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2915 } 2916 PetscCall(ISRestoreIndices(is[i],&idx)); 2917 PetscCall(ISDestroy(&is[i])); 2918 2919 k = 0; 2920 for (j=0; j<ov; j++) { /* for each overlap */ 2921 n = isz; 2922 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2923 row = nidx[k]; 2924 start = ai[row]; 2925 end = ai[row+1]; 2926 for (l = start; l<end; l++) { 2927 val = aj[l]; 2928 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2929 } 2930 } 2931 } 2932 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i))); 2933 } 2934 PetscCall(PetscBTDestroy(&table)); 2935 PetscCall(PetscFree(nidx)); 2936 PetscFunctionReturn(0); 2937 } 2938 2939 /* -------------------------------------------------------------- */ 2940 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2941 { 2942 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2943 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2944 const PetscInt *row,*col; 2945 PetscInt *cnew,j,*lens; 2946 IS icolp,irowp; 2947 PetscInt *cwork = NULL; 2948 PetscScalar *vwork = NULL; 2949 2950 PetscFunctionBegin; 2951 PetscCall(ISInvertPermutation(rowp,PETSC_DECIDE,&irowp)); 2952 PetscCall(ISGetIndices(irowp,&row)); 2953 PetscCall(ISInvertPermutation(colp,PETSC_DECIDE,&icolp)); 2954 PetscCall(ISGetIndices(icolp,&col)); 2955 2956 /* determine lengths of permuted rows */ 2957 PetscCall(PetscMalloc1(m+1,&lens)); 2958 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2959 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 2960 PetscCall(MatSetSizes(*B,m,n,m,n)); 2961 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 2962 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 2963 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens)); 2964 PetscCall(PetscFree(lens)); 2965 2966 PetscCall(PetscMalloc1(n,&cnew)); 2967 for (i=0; i<m; i++) { 2968 PetscCall(MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2969 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2970 PetscCall(MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES)); 2971 PetscCall(MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2972 } 2973 PetscCall(PetscFree(cnew)); 2974 2975 (*B)->assembled = PETSC_FALSE; 2976 2977 #if defined(PETSC_HAVE_DEVICE) 2978 PetscCall(MatBindToCPU(*B,A->boundtocpu)); 2979 #endif 2980 PetscCall(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 2981 PetscCall(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY)); 2982 PetscCall(ISRestoreIndices(irowp,&row)); 2983 PetscCall(ISRestoreIndices(icolp,&col)); 2984 PetscCall(ISDestroy(&irowp)); 2985 PetscCall(ISDestroy(&icolp)); 2986 if (rowp == colp) { 2987 PetscCall(MatPropagateSymmetryOptions(A,*B)); 2988 } 2989 PetscFunctionReturn(0); 2990 } 2991 2992 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2993 { 2994 PetscFunctionBegin; 2995 /* If the two matrices have the same copy implementation, use fast copy. */ 2996 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2997 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2998 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2999 const PetscScalar *aa; 3000 3001 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 3002 PetscCheck(a->i[A->rmap->n] == b->i[B->rmap->n],PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different %" PetscInt_FMT " != %" PetscInt_FMT,a->i[A->rmap->n],b->i[B->rmap->n]); 3003 PetscCall(PetscArraycpy(b->a,aa,a->i[A->rmap->n])); 3004 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 3005 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 3006 } else { 3007 PetscCall(MatCopy_Basic(A,B,str)); 3008 } 3009 PetscFunctionReturn(0); 3010 } 3011 3012 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 3013 { 3014 PetscFunctionBegin; 3015 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,NULL)); 3016 PetscFunctionReturn(0); 3017 } 3018 3019 PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 3020 { 3021 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3022 3023 PetscFunctionBegin; 3024 *array = a->a; 3025 PetscFunctionReturn(0); 3026 } 3027 3028 PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 3029 { 3030 PetscFunctionBegin; 3031 *array = NULL; 3032 PetscFunctionReturn(0); 3033 } 3034 3035 /* 3036 Computes the number of nonzeros per row needed for preallocation when X and Y 3037 have different nonzero structure. 3038 */ 3039 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 3040 { 3041 PetscInt i,j,k,nzx,nzy; 3042 3043 PetscFunctionBegin; 3044 /* Set the number of nonzeros in the new matrix */ 3045 for (i=0; i<m; i++) { 3046 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 3047 nzx = xi[i+1] - xi[i]; 3048 nzy = yi[i+1] - yi[i]; 3049 nnz[i] = 0; 3050 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 3051 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 3052 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 3053 nnz[i]++; 3054 } 3055 for (; k<nzy; k++) nnz[i]++; 3056 } 3057 PetscFunctionReturn(0); 3058 } 3059 3060 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 3061 { 3062 PetscInt m = Y->rmap->N; 3063 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 3064 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 3065 3066 PetscFunctionBegin; 3067 /* Set the number of nonzeros in the new matrix */ 3068 PetscCall(MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz)); 3069 PetscFunctionReturn(0); 3070 } 3071 3072 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 3073 { 3074 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 3075 3076 PetscFunctionBegin; 3077 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 3078 PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE; 3079 if (e) { 3080 PetscCall(PetscArraycmp(x->i,y->i,Y->rmap->n+1,&e)); 3081 if (e) { 3082 PetscCall(PetscArraycmp(x->j,y->j,y->nz,&e)); 3083 if (e) str = SAME_NONZERO_PATTERN; 3084 } 3085 } 3086 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatStructure is not SAME_NONZERO_PATTERN"); 3087 } 3088 if (str == SAME_NONZERO_PATTERN) { 3089 const PetscScalar *xa; 3090 PetscScalar *ya,alpha = a; 3091 PetscBLASInt one = 1,bnz; 3092 3093 PetscCall(PetscBLASIntCast(x->nz,&bnz)); 3094 PetscCall(MatSeqAIJGetArray(Y,&ya)); 3095 PetscCall(MatSeqAIJGetArrayRead(X,&xa)); 3096 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,xa,&one,ya,&one)); 3097 PetscCall(MatSeqAIJRestoreArrayRead(X,&xa)); 3098 PetscCall(MatSeqAIJRestoreArray(Y,&ya)); 3099 PetscCall(PetscLogFlops(2.0*bnz)); 3100 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 3101 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3102 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 3103 PetscCall(MatAXPY_Basic(Y,a,X,str)); 3104 } else { 3105 Mat B; 3106 PetscInt *nnz; 3107 PetscCall(PetscMalloc1(Y->rmap->N,&nnz)); 3108 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y),&B)); 3109 PetscCall(PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name)); 3110 PetscCall(MatSetLayouts(B,Y->rmap,Y->cmap)); 3111 PetscCall(MatSetType(B,((PetscObject)Y)->type_name)); 3112 PetscCall(MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz)); 3113 PetscCall(MatSeqAIJSetPreallocation(B,0,nnz)); 3114 PetscCall(MatAXPY_BasicWithPreallocation(B,Y,a,X,str)); 3115 PetscCall(MatHeaderMerge(Y,&B)); 3116 PetscCall(PetscFree(nnz)); 3117 } 3118 PetscFunctionReturn(0); 3119 } 3120 3121 PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 3122 { 3123 #if defined(PETSC_USE_COMPLEX) 3124 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3125 PetscInt i,nz; 3126 PetscScalar *a; 3127 3128 PetscFunctionBegin; 3129 nz = aij->nz; 3130 PetscCall(MatSeqAIJGetArray(mat,&a)); 3131 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 3132 PetscCall(MatSeqAIJRestoreArray(mat,&a)); 3133 #else 3134 PetscFunctionBegin; 3135 #endif 3136 PetscFunctionReturn(0); 3137 } 3138 3139 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3140 { 3141 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3142 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3143 PetscReal atmp; 3144 PetscScalar *x; 3145 const MatScalar *aa,*av; 3146 3147 PetscFunctionBegin; 3148 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3149 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3150 aa = av; 3151 ai = a->i; 3152 aj = a->j; 3153 3154 PetscCall(VecSet(v,0.0)); 3155 PetscCall(VecGetArrayWrite(v,&x)); 3156 PetscCall(VecGetLocalSize(v,&n)); 3157 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3158 for (i=0; i<m; i++) { 3159 ncols = ai[1] - ai[0]; ai++; 3160 for (j=0; j<ncols; j++) { 3161 atmp = PetscAbsScalar(*aa); 3162 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3163 aa++; aj++; 3164 } 3165 } 3166 PetscCall(VecRestoreArrayWrite(v,&x)); 3167 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3168 PetscFunctionReturn(0); 3169 } 3170 3171 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3172 { 3173 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3174 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3175 PetscScalar *x; 3176 const MatScalar *aa,*av; 3177 3178 PetscFunctionBegin; 3179 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3180 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3181 aa = av; 3182 ai = a->i; 3183 aj = a->j; 3184 3185 PetscCall(VecSet(v,0.0)); 3186 PetscCall(VecGetArrayWrite(v,&x)); 3187 PetscCall(VecGetLocalSize(v,&n)); 3188 PetscCheck(n == A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3189 for (i=0; i<m; i++) { 3190 ncols = ai[1] - ai[0]; ai++; 3191 if (ncols == A->cmap->n) { /* row is dense */ 3192 x[i] = *aa; if (idx) idx[i] = 0; 3193 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3194 x[i] = 0.0; 3195 if (idx) { 3196 for (j=0; j<ncols; j++) { /* find first implicit 0.0 in the row */ 3197 if (aj[j] > j) { 3198 idx[i] = j; 3199 break; 3200 } 3201 } 3202 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3203 if (j==ncols && j < A->cmap->n) idx[i] = j; 3204 } 3205 } 3206 for (j=0; j<ncols; j++) { 3207 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3208 aa++; aj++; 3209 } 3210 } 3211 PetscCall(VecRestoreArrayWrite(v,&x)); 3212 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3213 PetscFunctionReturn(0); 3214 } 3215 3216 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3217 { 3218 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3219 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3220 PetscScalar *x; 3221 const MatScalar *aa,*av; 3222 3223 PetscFunctionBegin; 3224 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3225 aa = av; 3226 ai = a->i; 3227 aj = a->j; 3228 3229 PetscCall(VecSet(v,0.0)); 3230 PetscCall(VecGetArrayWrite(v,&x)); 3231 PetscCall(VecGetLocalSize(v,&n)); 3232 PetscCheck(n == m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %" PetscInt_FMT " vs. %" PetscInt_FMT " rows", m, n); 3233 for (i=0; i<m; i++) { 3234 ncols = ai[1] - ai[0]; ai++; 3235 if (ncols == A->cmap->n) { /* row is dense */ 3236 x[i] = *aa; if (idx) idx[i] = 0; 3237 } else { /* row is sparse so already KNOW minimum is 0.0 or higher */ 3238 x[i] = 0.0; 3239 if (idx) { /* find first implicit 0.0 in the row */ 3240 for (j=0; j<ncols; j++) { 3241 if (aj[j] > j) { 3242 idx[i] = j; 3243 break; 3244 } 3245 } 3246 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3247 if (j==ncols && j < A->cmap->n) idx[i] = j; 3248 } 3249 } 3250 for (j=0; j<ncols; j++) { 3251 if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3252 aa++; aj++; 3253 } 3254 } 3255 PetscCall(VecRestoreArrayWrite(v,&x)); 3256 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3257 PetscFunctionReturn(0); 3258 } 3259 3260 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3261 { 3262 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3263 PetscInt i,j,m = A->rmap->n,ncols,n; 3264 const PetscInt *ai,*aj; 3265 PetscScalar *x; 3266 const MatScalar *aa,*av; 3267 3268 PetscFunctionBegin; 3269 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3270 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3271 aa = av; 3272 ai = a->i; 3273 aj = a->j; 3274 3275 PetscCall(VecSet(v,0.0)); 3276 PetscCall(VecGetArrayWrite(v,&x)); 3277 PetscCall(VecGetLocalSize(v,&n)); 3278 PetscCheck(n == m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3279 for (i=0; i<m; i++) { 3280 ncols = ai[1] - ai[0]; ai++; 3281 if (ncols == A->cmap->n) { /* row is dense */ 3282 x[i] = *aa; if (idx) idx[i] = 0; 3283 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3284 x[i] = 0.0; 3285 if (idx) { /* find first implicit 0.0 in the row */ 3286 for (j=0; j<ncols; j++) { 3287 if (aj[j] > j) { 3288 idx[i] = j; 3289 break; 3290 } 3291 } 3292 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3293 if (j==ncols && j < A->cmap->n) idx[i] = j; 3294 } 3295 } 3296 for (j=0; j<ncols; j++) { 3297 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3298 aa++; aj++; 3299 } 3300 } 3301 PetscCall(VecRestoreArrayWrite(v,&x)); 3302 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3303 PetscFunctionReturn(0); 3304 } 3305 3306 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3307 { 3308 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3309 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3310 MatScalar *diag,work[25],*v_work; 3311 const PetscReal shift = 0.0; 3312 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3313 3314 PetscFunctionBegin; 3315 allowzeropivot = PetscNot(A->erroriffailure); 3316 if (a->ibdiagvalid) { 3317 if (values) *values = a->ibdiag; 3318 PetscFunctionReturn(0); 3319 } 3320 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 3321 if (!a->ibdiag) { 3322 PetscCall(PetscMalloc1(bs2*mbs,&a->ibdiag)); 3323 PetscCall(PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar))); 3324 } 3325 diag = a->ibdiag; 3326 if (values) *values = a->ibdiag; 3327 /* factor and invert each block */ 3328 switch (bs) { 3329 case 1: 3330 for (i=0; i<mbs; i++) { 3331 PetscCall(MatGetValues(A,1,&i,1,&i,diag+i)); 3332 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3333 if (allowzeropivot) { 3334 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3335 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3336 A->factorerror_zeropivot_row = i; 3337 PetscCall(PetscInfo(A,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON)); 3338 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON); 3339 } 3340 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3341 } 3342 break; 3343 case 2: 3344 for (i=0; i<mbs; i++) { 3345 ij[0] = 2*i; ij[1] = 2*i + 1; 3346 PetscCall(MatGetValues(A,2,ij,2,ij,diag)); 3347 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 3348 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3349 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 3350 diag += 4; 3351 } 3352 break; 3353 case 3: 3354 for (i=0; i<mbs; i++) { 3355 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3356 PetscCall(MatGetValues(A,3,ij,3,ij,diag)); 3357 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 3358 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3359 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 3360 diag += 9; 3361 } 3362 break; 3363 case 4: 3364 for (i=0; i<mbs; i++) { 3365 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3366 PetscCall(MatGetValues(A,4,ij,4,ij,diag)); 3367 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 3368 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3369 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 3370 diag += 16; 3371 } 3372 break; 3373 case 5: 3374 for (i=0; i<mbs; i++) { 3375 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3376 PetscCall(MatGetValues(A,5,ij,5,ij,diag)); 3377 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 3378 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3379 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 3380 diag += 25; 3381 } 3382 break; 3383 case 6: 3384 for (i=0; i<mbs; i++) { 3385 ij[0] = 6*i; ij[1] = 6*i + 1; ij[2] = 6*i + 2; ij[3] = 6*i + 3; ij[4] = 6*i + 4; ij[5] = 6*i + 5; 3386 PetscCall(MatGetValues(A,6,ij,6,ij,diag)); 3387 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 3388 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3389 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 3390 diag += 36; 3391 } 3392 break; 3393 case 7: 3394 for (i=0; i<mbs; i++) { 3395 ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6; 3396 PetscCall(MatGetValues(A,7,ij,7,ij,diag)); 3397 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 3398 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3399 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 3400 diag += 49; 3401 } 3402 break; 3403 default: 3404 PetscCall(PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ)); 3405 for (i=0; i<mbs; i++) { 3406 for (j=0; j<bs; j++) { 3407 IJ[j] = bs*i + j; 3408 } 3409 PetscCall(MatGetValues(A,bs,IJ,bs,IJ,diag)); 3410 PetscCall(PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 3411 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3412 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bs)); 3413 diag += bs2; 3414 } 3415 PetscCall(PetscFree3(v_work,v_pivots,IJ)); 3416 } 3417 a->ibdiagvalid = PETSC_TRUE; 3418 PetscFunctionReturn(0); 3419 } 3420 3421 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3422 { 3423 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3424 PetscScalar a,*aa; 3425 PetscInt m,n,i,j,col; 3426 3427 PetscFunctionBegin; 3428 if (!x->assembled) { 3429 PetscCall(MatGetSize(x,&m,&n)); 3430 for (i=0; i<m; i++) { 3431 for (j=0; j<aij->imax[i]; j++) { 3432 PetscCall(PetscRandomGetValue(rctx,&a)); 3433 col = (PetscInt)(n*PetscRealPart(a)); 3434 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3435 } 3436 } 3437 } else { 3438 PetscCall(MatSeqAIJGetArrayWrite(x,&aa)); 3439 for (i=0; i<aij->nz; i++) PetscCall(PetscRandomGetValue(rctx,aa+i)); 3440 PetscCall(MatSeqAIJRestoreArrayWrite(x,&aa)); 3441 } 3442 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3443 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3444 PetscFunctionReturn(0); 3445 } 3446 3447 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3448 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx) 3449 { 3450 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3451 PetscScalar a; 3452 PetscInt m,n,i,j,col,nskip; 3453 3454 PetscFunctionBegin; 3455 nskip = high - low; 3456 PetscCall(MatGetSize(x,&m,&n)); 3457 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3458 for (i=0; i<m; i++) { 3459 for (j=0; j<aij->imax[i]; j++) { 3460 PetscCall(PetscRandomGetValue(rctx,&a)); 3461 col = (PetscInt)(n*PetscRealPart(a)); 3462 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3463 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3464 } 3465 } 3466 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3467 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3468 PetscFunctionReturn(0); 3469 } 3470 3471 /* -------------------------------------------------------------------*/ 3472 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3473 MatGetRow_SeqAIJ, 3474 MatRestoreRow_SeqAIJ, 3475 MatMult_SeqAIJ, 3476 /* 4*/ MatMultAdd_SeqAIJ, 3477 MatMultTranspose_SeqAIJ, 3478 MatMultTransposeAdd_SeqAIJ, 3479 NULL, 3480 NULL, 3481 NULL, 3482 /* 10*/ NULL, 3483 MatLUFactor_SeqAIJ, 3484 NULL, 3485 MatSOR_SeqAIJ, 3486 MatTranspose_SeqAIJ, 3487 /*1 5*/ MatGetInfo_SeqAIJ, 3488 MatEqual_SeqAIJ, 3489 MatGetDiagonal_SeqAIJ, 3490 MatDiagonalScale_SeqAIJ, 3491 MatNorm_SeqAIJ, 3492 /* 20*/ NULL, 3493 MatAssemblyEnd_SeqAIJ, 3494 MatSetOption_SeqAIJ, 3495 MatZeroEntries_SeqAIJ, 3496 /* 24*/ MatZeroRows_SeqAIJ, 3497 NULL, 3498 NULL, 3499 NULL, 3500 NULL, 3501 /* 29*/ MatSetUp_SeqAIJ, 3502 NULL, 3503 NULL, 3504 NULL, 3505 NULL, 3506 /* 34*/ MatDuplicate_SeqAIJ, 3507 NULL, 3508 NULL, 3509 MatILUFactor_SeqAIJ, 3510 NULL, 3511 /* 39*/ MatAXPY_SeqAIJ, 3512 MatCreateSubMatrices_SeqAIJ, 3513 MatIncreaseOverlap_SeqAIJ, 3514 MatGetValues_SeqAIJ, 3515 MatCopy_SeqAIJ, 3516 /* 44*/ MatGetRowMax_SeqAIJ, 3517 MatScale_SeqAIJ, 3518 MatShift_SeqAIJ, 3519 MatDiagonalSet_SeqAIJ, 3520 MatZeroRowsColumns_SeqAIJ, 3521 /* 49*/ MatSetRandom_SeqAIJ, 3522 MatGetRowIJ_SeqAIJ, 3523 MatRestoreRowIJ_SeqAIJ, 3524 MatGetColumnIJ_SeqAIJ, 3525 MatRestoreColumnIJ_SeqAIJ, 3526 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3527 NULL, 3528 NULL, 3529 MatPermute_SeqAIJ, 3530 NULL, 3531 /* 59*/ NULL, 3532 MatDestroy_SeqAIJ, 3533 MatView_SeqAIJ, 3534 NULL, 3535 NULL, 3536 /* 64*/ NULL, 3537 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3538 NULL, 3539 NULL, 3540 NULL, 3541 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3542 MatGetRowMinAbs_SeqAIJ, 3543 NULL, 3544 NULL, 3545 NULL, 3546 /* 74*/ NULL, 3547 MatFDColoringApply_AIJ, 3548 NULL, 3549 NULL, 3550 NULL, 3551 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3552 NULL, 3553 NULL, 3554 NULL, 3555 MatLoad_SeqAIJ, 3556 /* 84*/ MatIsSymmetric_SeqAIJ, 3557 MatIsHermitian_SeqAIJ, 3558 NULL, 3559 NULL, 3560 NULL, 3561 /* 89*/ NULL, 3562 NULL, 3563 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3564 NULL, 3565 NULL, 3566 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3567 NULL, 3568 NULL, 3569 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3570 NULL, 3571 /* 99*/ MatProductSetFromOptions_SeqAIJ, 3572 NULL, 3573 NULL, 3574 MatConjugate_SeqAIJ, 3575 NULL, 3576 /*104*/ MatSetValuesRow_SeqAIJ, 3577 MatRealPart_SeqAIJ, 3578 MatImaginaryPart_SeqAIJ, 3579 NULL, 3580 NULL, 3581 /*109*/ MatMatSolve_SeqAIJ, 3582 NULL, 3583 MatGetRowMin_SeqAIJ, 3584 NULL, 3585 MatMissingDiagonal_SeqAIJ, 3586 /*114*/ NULL, 3587 NULL, 3588 NULL, 3589 NULL, 3590 NULL, 3591 /*119*/ NULL, 3592 NULL, 3593 NULL, 3594 NULL, 3595 MatGetMultiProcBlock_SeqAIJ, 3596 /*124*/ MatFindNonzeroRows_SeqAIJ, 3597 MatGetColumnReductions_SeqAIJ, 3598 MatInvertBlockDiagonal_SeqAIJ, 3599 MatInvertVariableBlockDiagonal_SeqAIJ, 3600 NULL, 3601 /*129*/ NULL, 3602 NULL, 3603 NULL, 3604 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3605 MatTransposeColoringCreate_SeqAIJ, 3606 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3607 MatTransColoringApplyDenToSp_SeqAIJ, 3608 NULL, 3609 NULL, 3610 MatRARtNumeric_SeqAIJ_SeqAIJ, 3611 /*139*/NULL, 3612 NULL, 3613 NULL, 3614 MatFDColoringSetUp_SeqXAIJ, 3615 MatFindOffBlockDiagonalEntries_SeqAIJ, 3616 MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3617 /*145*/MatDestroySubMatrices_SeqAIJ, 3618 NULL, 3619 NULL 3620 }; 3621 3622 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3623 { 3624 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3625 PetscInt i,nz,n; 3626 3627 PetscFunctionBegin; 3628 nz = aij->maxnz; 3629 n = mat->rmap->n; 3630 for (i=0; i<nz; i++) { 3631 aij->j[i] = indices[i]; 3632 } 3633 aij->nz = nz; 3634 for (i=0; i<n; i++) { 3635 aij->ilen[i] = aij->imax[i]; 3636 } 3637 PetscFunctionReturn(0); 3638 } 3639 3640 /* 3641 * Given a sparse matrix with global column indices, compact it by using a local column space. 3642 * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3643 */ 3644 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3645 { 3646 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3647 PetscTable gid1_lid1; 3648 PetscTablePosition tpos; 3649 PetscInt gid,lid,i,ec,nz = aij->nz; 3650 PetscInt *garray,*jj = aij->j; 3651 3652 PetscFunctionBegin; 3653 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3654 PetscValidPointer(mapping,2); 3655 /* use a table */ 3656 PetscCall(PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1)); 3657 ec = 0; 3658 for (i=0; i<nz; i++) { 3659 PetscInt data,gid1 = jj[i] + 1; 3660 PetscCall(PetscTableFind(gid1_lid1,gid1,&data)); 3661 if (!data) { 3662 /* one based table */ 3663 PetscCall(PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES)); 3664 } 3665 } 3666 /* form array of columns we need */ 3667 PetscCall(PetscMalloc1(ec,&garray)); 3668 PetscCall(PetscTableGetHeadPosition(gid1_lid1,&tpos)); 3669 while (tpos) { 3670 PetscCall(PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid)); 3671 gid--; 3672 lid--; 3673 garray[lid] = gid; 3674 } 3675 PetscCall(PetscSortInt(ec,garray)); /* sort, and rebuild */ 3676 PetscCall(PetscTableRemoveAll(gid1_lid1)); 3677 for (i=0; i<ec; i++) { 3678 PetscCall(PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES)); 3679 } 3680 /* compact out the extra columns in B */ 3681 for (i=0; i<nz; i++) { 3682 PetscInt gid1 = jj[i] + 1; 3683 PetscCall(PetscTableFind(gid1_lid1,gid1,&lid)); 3684 lid--; 3685 jj[i] = lid; 3686 } 3687 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3688 PetscCall(PetscTableDestroy(&gid1_lid1)); 3689 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat),ec,ec,1,&mat->cmap)); 3690 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping)); 3691 PetscCall(ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH)); 3692 PetscFunctionReturn(0); 3693 } 3694 3695 /*@ 3696 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3697 in the matrix. 3698 3699 Input Parameters: 3700 + mat - the SeqAIJ matrix 3701 - indices - the column indices 3702 3703 Level: advanced 3704 3705 Notes: 3706 This can be called if you have precomputed the nonzero structure of the 3707 matrix and want to provide it to the matrix object to improve the performance 3708 of the MatSetValues() operation. 3709 3710 You MUST have set the correct numbers of nonzeros per row in the call to 3711 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3712 3713 MUST be called before any calls to MatSetValues(); 3714 3715 The indices should start with zero, not one. 3716 3717 @*/ 3718 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3719 { 3720 PetscFunctionBegin; 3721 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3722 PetscValidIntPointer(indices,2); 3723 PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices)); 3724 PetscFunctionReturn(0); 3725 } 3726 3727 /* ----------------------------------------------------------------------------------------*/ 3728 3729 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3730 { 3731 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3732 size_t nz = aij->i[mat->rmap->n]; 3733 3734 PetscFunctionBegin; 3735 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3736 3737 /* allocate space for values if not already there */ 3738 if (!aij->saved_values) { 3739 PetscCall(PetscMalloc1(nz+1,&aij->saved_values)); 3740 PetscCall(PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar))); 3741 } 3742 3743 /* copy values over */ 3744 PetscCall(PetscArraycpy(aij->saved_values,aij->a,nz)); 3745 PetscFunctionReturn(0); 3746 } 3747 3748 /*@ 3749 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3750 example, reuse of the linear part of a Jacobian, while recomputing the 3751 nonlinear portion. 3752 3753 Collect on Mat 3754 3755 Input Parameters: 3756 . mat - the matrix (currently only AIJ matrices support this option) 3757 3758 Level: advanced 3759 3760 Common Usage, with SNESSolve(): 3761 $ Create Jacobian matrix 3762 $ Set linear terms into matrix 3763 $ Apply boundary conditions to matrix, at this time matrix must have 3764 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3765 $ boundary conditions again will not change the nonzero structure 3766 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3767 $ ierr = MatStoreValues(mat); 3768 $ Call SNESSetJacobian() with matrix 3769 $ In your Jacobian routine 3770 $ ierr = MatRetrieveValues(mat); 3771 $ Set nonlinear terms in matrix 3772 3773 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3774 $ // build linear portion of Jacobian 3775 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3776 $ ierr = MatStoreValues(mat); 3777 $ loop over nonlinear iterations 3778 $ ierr = MatRetrieveValues(mat); 3779 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3780 $ // call MatAssemblyBegin/End() on matrix 3781 $ Solve linear system with Jacobian 3782 $ endloop 3783 3784 Notes: 3785 Matrix must already be assemblied before calling this routine 3786 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3787 calling this routine. 3788 3789 When this is called multiple times it overwrites the previous set of stored values 3790 and does not allocated additional space. 3791 3792 .seealso: `MatRetrieveValues()` 3793 3794 @*/ 3795 PetscErrorCode MatStoreValues(Mat mat) 3796 { 3797 PetscFunctionBegin; 3798 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3799 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3800 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3801 PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat)); 3802 PetscFunctionReturn(0); 3803 } 3804 3805 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3806 { 3807 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3808 PetscInt nz = aij->i[mat->rmap->n]; 3809 3810 PetscFunctionBegin; 3811 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3812 PetscCheck(aij->saved_values,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3813 /* copy values over */ 3814 PetscCall(PetscArraycpy(aij->a,aij->saved_values,nz)); 3815 PetscFunctionReturn(0); 3816 } 3817 3818 /*@ 3819 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3820 example, reuse of the linear part of a Jacobian, while recomputing the 3821 nonlinear portion. 3822 3823 Collect on Mat 3824 3825 Input Parameters: 3826 . mat - the matrix (currently only AIJ matrices support this option) 3827 3828 Level: advanced 3829 3830 .seealso: `MatStoreValues()` 3831 3832 @*/ 3833 PetscErrorCode MatRetrieveValues(Mat mat) 3834 { 3835 PetscFunctionBegin; 3836 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3837 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3838 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3839 PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat)); 3840 PetscFunctionReturn(0); 3841 } 3842 3843 /* --------------------------------------------------------------------------------*/ 3844 /*@C 3845 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3846 (the default parallel PETSc format). For good matrix assembly performance 3847 the user should preallocate the matrix storage by setting the parameter nz 3848 (or the array nnz). By setting these parameters accurately, performance 3849 during matrix assembly can be increased by more than a factor of 50. 3850 3851 Collective 3852 3853 Input Parameters: 3854 + comm - MPI communicator, set to PETSC_COMM_SELF 3855 . m - number of rows 3856 . n - number of columns 3857 . nz - number of nonzeros per row (same for all rows) 3858 - nnz - array containing the number of nonzeros in the various rows 3859 (possibly different for each row) or NULL 3860 3861 Output Parameter: 3862 . A - the matrix 3863 3864 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3865 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3866 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3867 3868 Notes: 3869 If nnz is given then nz is ignored 3870 3871 The AIJ format (also called the Yale sparse matrix format or 3872 compressed row storage), is fully compatible with standard Fortran 77 3873 storage. That is, the stored row and column indices can begin at 3874 either one (as in Fortran) or zero. See the users' manual for details. 3875 3876 Specify the preallocated storage with either nz or nnz (not both). 3877 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3878 allocation. For large problems you MUST preallocate memory or you 3879 will get TERRIBLE performance, see the users' manual chapter on matrices. 3880 3881 By default, this format uses inodes (identical nodes) when possible, to 3882 improve numerical efficiency of matrix-vector products and solves. We 3883 search for consecutive rows with the same nonzero structure, thereby 3884 reusing matrix information to achieve increased efficiency. 3885 3886 Options Database Keys: 3887 + -mat_no_inode - Do not use inodes 3888 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3889 3890 Level: intermediate 3891 3892 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()` 3893 3894 @*/ 3895 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3896 { 3897 PetscFunctionBegin; 3898 PetscCall(MatCreate(comm,A)); 3899 PetscCall(MatSetSizes(*A,m,n,m,n)); 3900 PetscCall(MatSetType(*A,MATSEQAIJ)); 3901 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz)); 3902 PetscFunctionReturn(0); 3903 } 3904 3905 /*@C 3906 MatSeqAIJSetPreallocation - For good matrix assembly performance 3907 the user should preallocate the matrix storage by setting the parameter nz 3908 (or the array nnz). By setting these parameters accurately, performance 3909 during matrix assembly can be increased by more than a factor of 50. 3910 3911 Collective 3912 3913 Input Parameters: 3914 + B - The matrix 3915 . nz - number of nonzeros per row (same for all rows) 3916 - nnz - array containing the number of nonzeros in the various rows 3917 (possibly different for each row) or NULL 3918 3919 Notes: 3920 If nnz is given then nz is ignored 3921 3922 The AIJ format (also called the Yale sparse matrix format or 3923 compressed row storage), is fully compatible with standard Fortran 77 3924 storage. That is, the stored row and column indices can begin at 3925 either one (as in Fortran) or zero. See the users' manual for details. 3926 3927 Specify the preallocated storage with either nz or nnz (not both). 3928 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3929 allocation. For large problems you MUST preallocate memory or you 3930 will get TERRIBLE performance, see the users' manual chapter on matrices. 3931 3932 You can call MatGetInfo() to get information on how effective the preallocation was; 3933 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3934 You can also run with the option -info and look for messages with the string 3935 malloc in them to see if additional memory allocation was needed. 3936 3937 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3938 entries or columns indices 3939 3940 By default, this format uses inodes (identical nodes) when possible, to 3941 improve numerical efficiency of matrix-vector products and solves. We 3942 search for consecutive rows with the same nonzero structure, thereby 3943 reusing matrix information to achieve increased efficiency. 3944 3945 Options Database Keys: 3946 + -mat_no_inode - Do not use inodes 3947 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3948 3949 Level: intermediate 3950 3951 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatSetValues()`, `MatSeqAIJSetColumnIndices()`, `MatCreateSeqAIJWithArrays()`, `MatGetInfo()`, 3952 `MatSeqAIJSetTotalPreallocation()` 3953 3954 @*/ 3955 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3956 { 3957 PetscFunctionBegin; 3958 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3959 PetscValidType(B,1); 3960 PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz)); 3961 PetscFunctionReturn(0); 3962 } 3963 3964 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3965 { 3966 Mat_SeqAIJ *b; 3967 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3968 PetscInt i; 3969 3970 PetscFunctionBegin; 3971 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3972 if (nz == MAT_SKIP_ALLOCATION) { 3973 skipallocation = PETSC_TRUE; 3974 nz = 0; 3975 } 3976 PetscCall(PetscLayoutSetUp(B->rmap)); 3977 PetscCall(PetscLayoutSetUp(B->cmap)); 3978 3979 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3980 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %" PetscInt_FMT,nz); 3981 if (PetscUnlikelyDebug(nnz)) { 3982 for (i=0; i<B->rmap->n; i++) { 3983 PetscCheck(nnz[i] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT,i,nnz[i]); 3984 PetscCheck(nnz[i] <= B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " rowlength %" PetscInt_FMT,i,nnz[i],B->cmap->n); 3985 } 3986 } 3987 3988 B->preallocated = PETSC_TRUE; 3989 3990 b = (Mat_SeqAIJ*)B->data; 3991 3992 if (!skipallocation) { 3993 if (!b->imax) { 3994 PetscCall(PetscMalloc1(B->rmap->n,&b->imax)); 3995 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 3996 } 3997 if (!b->ilen) { 3998 /* b->ilen will count nonzeros in each row so far. */ 3999 PetscCall(PetscCalloc1(B->rmap->n,&b->ilen)); 4000 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4001 } else { 4002 PetscCall(PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt))); 4003 } 4004 if (!b->ipre) { 4005 PetscCall(PetscMalloc1(B->rmap->n,&b->ipre)); 4006 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4007 } 4008 if (!nnz) { 4009 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 4010 else if (nz < 0) nz = 1; 4011 nz = PetscMin(nz,B->cmap->n); 4012 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 4013 nz = nz*B->rmap->n; 4014 } else { 4015 PetscInt64 nz64 = 0; 4016 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 4017 PetscCall(PetscIntCast(nz64,&nz)); 4018 } 4019 4020 /* allocate the matrix space */ 4021 /* FIXME: should B's old memory be unlogged? */ 4022 PetscCall(MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i)); 4023 if (B->structure_only) { 4024 PetscCall(PetscMalloc1(nz,&b->j)); 4025 PetscCall(PetscMalloc1(B->rmap->n+1,&b->i)); 4026 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt))); 4027 } else { 4028 PetscCall(PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i)); 4029 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)))); 4030 } 4031 b->i[0] = 0; 4032 for (i=1; i<B->rmap->n+1; i++) { 4033 b->i[i] = b->i[i-1] + b->imax[i-1]; 4034 } 4035 if (B->structure_only) { 4036 b->singlemalloc = PETSC_FALSE; 4037 b->free_a = PETSC_FALSE; 4038 } else { 4039 b->singlemalloc = PETSC_TRUE; 4040 b->free_a = PETSC_TRUE; 4041 } 4042 b->free_ij = PETSC_TRUE; 4043 } else { 4044 b->free_a = PETSC_FALSE; 4045 b->free_ij = PETSC_FALSE; 4046 } 4047 4048 if (b->ipre && nnz != b->ipre && b->imax) { 4049 /* reserve user-requested sparsity */ 4050 PetscCall(PetscArraycpy(b->ipre,b->imax,B->rmap->n)); 4051 } 4052 4053 b->nz = 0; 4054 b->maxnz = nz; 4055 B->info.nz_unneeded = (double)b->maxnz; 4056 if (realalloc) { 4057 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 4058 } 4059 B->was_assembled = PETSC_FALSE; 4060 B->assembled = PETSC_FALSE; 4061 /* We simply deem preallocation has changed nonzero state. Updating the state 4062 will give clients (like AIJKokkos) a chance to know something has happened. 4063 */ 4064 B->nonzerostate++; 4065 PetscFunctionReturn(0); 4066 } 4067 4068 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4069 { 4070 Mat_SeqAIJ *a; 4071 PetscInt i; 4072 4073 PetscFunctionBegin; 4074 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4075 4076 /* Check local size. If zero, then return */ 4077 if (!A->rmap->n) PetscFunctionReturn(0); 4078 4079 a = (Mat_SeqAIJ*)A->data; 4080 /* if no saved info, we error out */ 4081 PetscCheck(a->ipre,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"No saved preallocation info "); 4082 4083 PetscCheck(a->i && a->j && a->a && a->imax && a->ilen,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Memory info is incomplete, and can not reset preallocation "); 4084 4085 PetscCall(PetscArraycpy(a->imax,a->ipre,A->rmap->n)); 4086 PetscCall(PetscArrayzero(a->ilen,A->rmap->n)); 4087 a->i[0] = 0; 4088 for (i=1; i<A->rmap->n+1; i++) { 4089 a->i[i] = a->i[i-1] + a->imax[i-1]; 4090 } 4091 A->preallocated = PETSC_TRUE; 4092 a->nz = 0; 4093 a->maxnz = a->i[A->rmap->n]; 4094 A->info.nz_unneeded = (double)a->maxnz; 4095 A->was_assembled = PETSC_FALSE; 4096 A->assembled = PETSC_FALSE; 4097 PetscFunctionReturn(0); 4098 } 4099 4100 /*@ 4101 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 4102 4103 Input Parameters: 4104 + B - the matrix 4105 . i - the indices into j for the start of each row (starts with zero) 4106 . j - the column indices for each row (starts with zero) these must be sorted for each row 4107 - v - optional values in the matrix 4108 4109 Level: developer 4110 4111 Notes: 4112 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 4113 4114 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4115 structure will be the union of all the previous nonzero structures. 4116 4117 Developer Notes: 4118 An optimization could be added to the implementation where it checks if the i, and j are identical to the current i and j and 4119 then just copies the v values directly with PetscMemcpy(). 4120 4121 This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them. 4122 4123 .seealso: `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatSeqAIJSetPreallocation()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MatResetPreallocation()` 4124 @*/ 4125 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 4126 { 4127 PetscFunctionBegin; 4128 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 4129 PetscValidType(B,1); 4130 PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v)); 4131 PetscFunctionReturn(0); 4132 } 4133 4134 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 4135 { 4136 PetscInt i; 4137 PetscInt m,n; 4138 PetscInt nz; 4139 PetscInt *nnz; 4140 4141 PetscFunctionBegin; 4142 PetscCheck(Ii[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4143 4144 PetscCall(PetscLayoutSetUp(B->rmap)); 4145 PetscCall(PetscLayoutSetUp(B->cmap)); 4146 4147 PetscCall(MatGetSize(B, &m, &n)); 4148 PetscCall(PetscMalloc1(m+1, &nnz)); 4149 for (i = 0; i < m; i++) { 4150 nz = Ii[i+1]- Ii[i]; 4151 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4152 nnz[i] = nz; 4153 } 4154 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4155 PetscCall(PetscFree(nnz)); 4156 4157 for (i = 0; i < m; i++) { 4158 PetscCall(MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES)); 4159 } 4160 4161 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 4162 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 4163 4164 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 4165 PetscFunctionReturn(0); 4166 } 4167 4168 /*@ 4169 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4170 4171 Input Parameters: 4172 + A - left-hand side matrix 4173 . B - right-hand side matrix 4174 - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 4175 4176 Output Parameter: 4177 . C - Kronecker product of A and B 4178 4179 Level: intermediate 4180 4181 Notes: 4182 MAT_REUSE_MATRIX can only be used when the nonzero structure of the product matrix has not changed from that last call to MatSeqAIJKron(). 4183 4184 .seealso: `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATKAIJ`, `MatReuse` 4185 @*/ 4186 PetscErrorCode MatSeqAIJKron(Mat A,Mat B,MatReuse reuse,Mat *C) 4187 { 4188 PetscFunctionBegin; 4189 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4190 PetscValidType(A,1); 4191 PetscValidHeaderSpecific(B,MAT_CLASSID,2); 4192 PetscValidType(B,2); 4193 PetscValidPointer(C,4); 4194 if (reuse == MAT_REUSE_MATRIX) { 4195 PetscValidHeaderSpecific(*C,MAT_CLASSID,4); 4196 PetscValidType(*C,4); 4197 } 4198 PetscTryMethod(A,"MatSeqAIJKron_C",(Mat,Mat,MatReuse,Mat*),(A,B,reuse,C)); 4199 PetscFunctionReturn(0); 4200 } 4201 4202 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A,Mat B,MatReuse reuse,Mat *C) 4203 { 4204 Mat newmat; 4205 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4206 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 4207 PetscScalar *v; 4208 const PetscScalar *aa,*ba; 4209 PetscInt *i,*j,m,n,p,q,nnz = 0,am = A->rmap->n,bm = B->rmap->n,an = A->cmap->n, bn = B->cmap->n; 4210 PetscBool flg; 4211 4212 PetscFunctionBegin; 4213 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4214 PetscCheck(A->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4215 PetscCheck(!B->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4216 PetscCheck(B->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4217 PetscCall(PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&flg)); 4218 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatType %s",((PetscObject)B)->type_name); 4219 PetscCheck(reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatReuse %d",(int)reuse); 4220 if (reuse == MAT_INITIAL_MATRIX) { 4221 PetscCall(PetscMalloc2(am*bm+1,&i,a->i[am]*b->i[bm],&j)); 4222 PetscCall(MatCreate(PETSC_COMM_SELF,&newmat)); 4223 PetscCall(MatSetSizes(newmat,am*bm,an*bn,am*bm,an*bn)); 4224 PetscCall(MatSetType(newmat,MATAIJ)); 4225 i[0] = 0; 4226 for (m = 0; m < am; ++m) { 4227 for (p = 0; p < bm; ++p) { 4228 i[m*bm + p + 1] = i[m*bm + p] + (a->i[m+1] - a->i[m]) * (b->i[p+1] - b->i[p]); 4229 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4230 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4231 j[nnz++] = a->j[n]*bn + b->j[q]; 4232 } 4233 } 4234 } 4235 } 4236 PetscCall(MatSeqAIJSetPreallocationCSR(newmat,i,j,NULL)); 4237 *C = newmat; 4238 PetscCall(PetscFree2(i,j)); 4239 nnz = 0; 4240 } 4241 PetscCall(MatSeqAIJGetArray(*C,&v)); 4242 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4243 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 4244 for (m = 0; m < am; ++m) { 4245 for (p = 0; p < bm; ++p) { 4246 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4247 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4248 v[nnz++] = aa[n] * ba[q]; 4249 } 4250 } 4251 } 4252 } 4253 PetscCall(MatSeqAIJRestoreArray(*C,&v)); 4254 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 4255 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 4256 PetscFunctionReturn(0); 4257 } 4258 4259 #include <../src/mat/impls/dense/seq/dense.h> 4260 #include <petsc/private/kernels/petscaxpy.h> 4261 4262 /* 4263 Computes (B'*A')' since computing B*A directly is untenable 4264 4265 n p p 4266 [ ] [ ] [ ] 4267 m [ A ] * n [ B ] = m [ C ] 4268 [ ] [ ] [ ] 4269 4270 */ 4271 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 4272 { 4273 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 4274 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 4275 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 4276 PetscInt i,j,n,m,q,p; 4277 const PetscInt *ii,*idx; 4278 const PetscScalar *b,*a,*a_q; 4279 PetscScalar *c,*c_q; 4280 PetscInt clda = sub_c->lda; 4281 PetscInt alda = sub_a->lda; 4282 4283 PetscFunctionBegin; 4284 m = A->rmap->n; 4285 n = A->cmap->n; 4286 p = B->cmap->n; 4287 a = sub_a->v; 4288 b = sub_b->a; 4289 c = sub_c->v; 4290 if (clda == m) { 4291 PetscCall(PetscArrayzero(c,m*p)); 4292 } else { 4293 for (j=0;j<p;j++) 4294 for (i=0;i<m;i++) 4295 c[j*clda + i] = 0.0; 4296 } 4297 ii = sub_b->i; 4298 idx = sub_b->j; 4299 for (i=0; i<n; i++) { 4300 q = ii[i+1] - ii[i]; 4301 while (q-->0) { 4302 c_q = c + clda*(*idx); 4303 a_q = a + alda*i; 4304 PetscKernelAXPY(c_q,*b,a_q,m); 4305 idx++; 4306 b++; 4307 } 4308 } 4309 PetscFunctionReturn(0); 4310 } 4311 4312 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C) 4313 { 4314 PetscInt m=A->rmap->n,n=B->cmap->n; 4315 PetscBool cisdense; 4316 4317 PetscFunctionBegin; 4318 PetscCheck(A->cmap->n == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %" PetscInt_FMT " != B->rmap->n %" PetscInt_FMT,A->cmap->n,B->rmap->n); 4319 PetscCall(MatSetSizes(C,m,n,m,n)); 4320 PetscCall(MatSetBlockSizesFromMats(C,A,B)); 4321 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 4322 if (!cisdense) { 4323 PetscCall(MatSetType(C,MATDENSE)); 4324 } 4325 PetscCall(MatSetUp(C)); 4326 4327 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4328 PetscFunctionReturn(0); 4329 } 4330 4331 /* ----------------------------------------------------------------*/ 4332 /*MC 4333 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4334 based on compressed sparse row format. 4335 4336 Options Database Keys: 4337 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4338 4339 Level: beginner 4340 4341 Notes: 4342 MatSetValues() may be called for this matrix type with a NULL argument for the numerical values, 4343 in this case the values associated with the rows and columns one passes in are set to zero 4344 in the matrix 4345 4346 MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no 4347 space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored 4348 4349 Developer Notes: 4350 It would be nice if all matrix formats supported passing NULL in for the numerical values 4351 4352 .seealso: `MatCreateSeqAIJ()`, `MatSetFromOptions()`, `MatSetType()`, `MatCreate()`, `MatType`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4353 M*/ 4354 4355 /*MC 4356 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4357 4358 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4359 and MATMPIAIJ otherwise. As a result, for single process communicators, 4360 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4361 for communicators controlling multiple processes. It is recommended that you call both of 4362 the above preallocation routines for simplicity. 4363 4364 Options Database Keys: 4365 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4366 4367 Developer Notes: 4368 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4369 enough exist. 4370 4371 Level: beginner 4372 4373 .seealso: `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MATSEQAIJ`, `MATMPIAIJ`, `MATSELL`, `MATSEQSELL`, `MATMPISELL` 4374 M*/ 4375 4376 /*MC 4377 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4378 4379 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4380 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4381 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4382 for communicators controlling multiple processes. It is recommended that you call both of 4383 the above preallocation routines for simplicity. 4384 4385 Options Database Keys: 4386 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4387 4388 Level: beginner 4389 4390 .seealso: `MatCreateMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL`, `MATSEQAIJCRL`, `MATMPIAIJCRL` 4391 M*/ 4392 4393 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4394 #if defined(PETSC_HAVE_ELEMENTAL) 4395 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4396 #endif 4397 #if defined(PETSC_HAVE_SCALAPACK) 4398 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*); 4399 #endif 4400 #if defined(PETSC_HAVE_HYPRE) 4401 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4402 #endif 4403 4404 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4405 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4406 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4407 4408 /*@C 4409 MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored 4410 4411 Not Collective 4412 4413 Input Parameter: 4414 . mat - a MATSEQAIJ matrix 4415 4416 Output Parameter: 4417 . array - pointer to the data 4418 4419 Level: intermediate 4420 4421 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4422 @*/ 4423 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4424 { 4425 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4426 4427 PetscFunctionBegin; 4428 if (aij->ops->getarray) { 4429 PetscCall((*aij->ops->getarray)(A,array)); 4430 } else { 4431 *array = aij->a; 4432 } 4433 PetscFunctionReturn(0); 4434 } 4435 4436 /*@C 4437 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4438 4439 Not Collective 4440 4441 Input Parameters: 4442 + mat - a MATSEQAIJ matrix 4443 - array - pointer to the data 4444 4445 Level: intermediate 4446 4447 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayF90()` 4448 @*/ 4449 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4450 { 4451 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4452 4453 PetscFunctionBegin; 4454 if (aij->ops->restorearray) { 4455 PetscCall((*aij->ops->restorearray)(A,array)); 4456 } else { 4457 *array = NULL; 4458 } 4459 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4460 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4461 PetscFunctionReturn(0); 4462 } 4463 4464 /*@C 4465 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored 4466 4467 Not Collective 4468 4469 Input Parameter: 4470 . mat - a MATSEQAIJ matrix 4471 4472 Output Parameter: 4473 . array - pointer to the data 4474 4475 Level: intermediate 4476 4477 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4478 @*/ 4479 PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array) 4480 { 4481 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4482 4483 PetscFunctionBegin; 4484 if (aij->ops->getarrayread) { 4485 PetscCall((*aij->ops->getarrayread)(A,array)); 4486 } else { 4487 *array = aij->a; 4488 } 4489 PetscFunctionReturn(0); 4490 } 4491 4492 /*@C 4493 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4494 4495 Not Collective 4496 4497 Input Parameter: 4498 . mat - a MATSEQAIJ matrix 4499 4500 Output Parameter: 4501 . array - pointer to the data 4502 4503 Level: intermediate 4504 4505 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4506 @*/ 4507 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array) 4508 { 4509 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4510 4511 PetscFunctionBegin; 4512 if (aij->ops->restorearrayread) { 4513 PetscCall((*aij->ops->restorearrayread)(A,array)); 4514 } else { 4515 *array = NULL; 4516 } 4517 PetscFunctionReturn(0); 4518 } 4519 4520 /*@C 4521 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a MATSEQAIJ matrix is stored 4522 4523 Not Collective 4524 4525 Input Parameter: 4526 . mat - a MATSEQAIJ matrix 4527 4528 Output Parameter: 4529 . array - pointer to the data 4530 4531 Level: intermediate 4532 4533 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArrayRead()` 4534 @*/ 4535 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A,PetscScalar **array) 4536 { 4537 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4538 4539 PetscFunctionBegin; 4540 if (aij->ops->getarraywrite) { 4541 PetscCall((*aij->ops->getarraywrite)(A,array)); 4542 } else { 4543 *array = aij->a; 4544 } 4545 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4546 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4547 PetscFunctionReturn(0); 4548 } 4549 4550 /*@C 4551 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4552 4553 Not Collective 4554 4555 Input Parameter: 4556 . mat - a MATSEQAIJ matrix 4557 4558 Output Parameter: 4559 . array - pointer to the data 4560 4561 Level: intermediate 4562 4563 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4564 @*/ 4565 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A,PetscScalar **array) 4566 { 4567 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4568 4569 PetscFunctionBegin; 4570 if (aij->ops->restorearraywrite) { 4571 PetscCall((*aij->ops->restorearraywrite)(A,array)); 4572 } else { 4573 *array = NULL; 4574 } 4575 PetscFunctionReturn(0); 4576 } 4577 4578 /*@C 4579 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the SEQAIJ matrix 4580 4581 Not Collective 4582 4583 Input Parameter: 4584 . mat - a matrix of type MATSEQAIJ or its subclasses 4585 4586 Output Parameters: 4587 + i - row map array of the matrix 4588 . j - column index array of the matrix 4589 . a - data array of the matrix 4590 - memtype - memory type of the arrays 4591 4592 Notes: 4593 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4594 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4595 4596 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4597 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4598 4599 Level: Developer 4600 4601 .seealso: `MatSeqAIJGetArray()`, `MatSeqAIJGetArrayRead()` 4602 @*/ 4603 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat,const PetscInt **i,const PetscInt **j,PetscScalar **a,PetscMemType *mtype) 4604 { 4605 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 4606 4607 PetscFunctionBegin; 4608 PetscCheck(mat->preallocated,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"matrix is not preallocated"); 4609 if (aij->ops->getcsrandmemtype) { 4610 PetscCall((*aij->ops->getcsrandmemtype)(mat,i,j,a,mtype)); 4611 } else { 4612 if (i) *i = aij->i; 4613 if (j) *j = aij->j; 4614 if (a) *a = aij->a; 4615 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4616 } 4617 PetscFunctionReturn(0); 4618 } 4619 4620 /*@C 4621 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4622 4623 Not Collective 4624 4625 Input Parameter: 4626 . mat - a MATSEQAIJ matrix 4627 4628 Output Parameter: 4629 . nz - the maximum number of nonzeros in any row 4630 4631 Level: intermediate 4632 4633 .seealso: `MatSeqAIJRestoreArray()`, `MatSeqAIJGetArrayF90()` 4634 @*/ 4635 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4636 { 4637 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4638 4639 PetscFunctionBegin; 4640 *nz = aij->rmax; 4641 PetscFunctionReturn(0); 4642 } 4643 4644 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, const PetscInt coo_i[], const PetscInt coo_j[]) 4645 { 4646 MPI_Comm comm; 4647 PetscInt *i,*j; 4648 PetscInt M,N,row; 4649 PetscCount k,p,q,nneg,nnz,start,end; /* Index the coo array, so use PetscCount as their type */ 4650 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4651 PetscInt *Aj; 4652 PetscScalar *Aa; 4653 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ*)(mat->data); 4654 MatType rtype; 4655 PetscCount *perm,*jmap; 4656 4657 PetscFunctionBegin; 4658 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4659 PetscCall(PetscObjectGetComm((PetscObject)mat,&comm)); 4660 PetscCall(MatGetSize(mat,&M,&N)); 4661 PetscCall(PetscMalloc2(coo_n,&i,coo_n,&j)); 4662 PetscCall(PetscArraycpy(i,coo_i,coo_n)); /* Make a copy since we'll modify it */ 4663 PetscCall(PetscArraycpy(j,coo_j,coo_n)); 4664 PetscCall(PetscMalloc1(coo_n,&perm)); 4665 for (k=0; k<coo_n; k++) { /* Ignore entries with negative row or col indices */ 4666 if (j[k] < 0) i[k] = -1; 4667 perm[k] = k; 4668 } 4669 4670 /* Sort by row */ 4671 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n,i,j,perm)); 4672 for (k=0; k<coo_n; k++) {if (i[k] >= 0) break;} /* Advance k to the first row with a non-negative index */ 4673 nneg = k; 4674 PetscCall(PetscMalloc1(coo_n-nneg+1,&jmap)); /* +1 to make a CSR-like data structure. jmap[i] originally is the number of repeats for i-th nonzero */ 4675 nnz = 0; /* Total number of unique nonzeros to be counted */ 4676 jmap++; /* Inc jmap by 1 for convinience */ 4677 4678 PetscCall(PetscCalloc1(M+1,&Ai)); /* CSR of A */ 4679 PetscCall(PetscMalloc1(coo_n-nneg,&Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4680 4681 /* In each row, sort by column, then unique column indices to get row length */ 4682 Ai++; /* Inc by 1 for convinience */ 4683 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4684 while (k<coo_n) { 4685 row = i[k]; 4686 start = k; /* [start,end) indices for this row */ 4687 while (k<coo_n && i[k] == row) k++; 4688 end = k; 4689 PetscCall(PetscSortIntWithCountArray(end-start,j+start,perm+start)); 4690 /* Find number of unique col entries in this row */ 4691 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4692 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4693 Ai[row] = 1; 4694 nnz++; 4695 4696 for (p=start+1; p<end; p++) { /* Scan remaining nonzero in this row */ 4697 if (j[p] != j[p-1]) { /* Meet a new nonzero */ 4698 q++; 4699 jmap[q] = 1; 4700 Aj[q] = j[p]; 4701 Ai[row]++; 4702 nnz++; 4703 } else { 4704 jmap[q]++; 4705 } 4706 } 4707 q++; /* Move to next row and thus next unique nonzero */ 4708 } 4709 PetscCall(PetscFree2(i,j)); 4710 4711 Ai--; /* Back to the beginning of Ai[] */ 4712 for (k=0; k<M; k++) Ai[k+1] += Ai[k]; 4713 jmap--; /* Back to the beginning of jmap[] */ 4714 jmap[0] = 0; 4715 for (k=0; k<nnz; k++) jmap[k+1] += jmap[k]; 4716 if (nnz < coo_n-nneg) { /* Realloc with actual number of unique nonzeros */ 4717 PetscCount *jmap_new; 4718 PetscInt *Aj_new; 4719 4720 PetscCall(PetscMalloc1(nnz+1,&jmap_new)); 4721 PetscCall(PetscArraycpy(jmap_new,jmap,nnz+1)); 4722 PetscCall(PetscFree(jmap)); 4723 jmap = jmap_new; 4724 4725 PetscCall(PetscMalloc1(nnz,&Aj_new)); 4726 PetscCall(PetscArraycpy(Aj_new,Aj,nnz)); 4727 PetscCall(PetscFree(Aj)); 4728 Aj = Aj_new; 4729 } 4730 4731 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4732 PetscCount *perm_new; 4733 4734 PetscCall(PetscMalloc1(coo_n-nneg,&perm_new)); 4735 PetscCall(PetscArraycpy(perm_new,perm+nneg,coo_n-nneg)); 4736 PetscCall(PetscFree(perm)); 4737 perm = perm_new; 4738 } 4739 4740 PetscCall(MatGetRootType_Private(mat,&rtype)); 4741 PetscCall(PetscCalloc1(nnz,&Aa)); /* Zero the matrix */ 4742 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF,M,N,Ai,Aj,Aa,rtype,mat)); 4743 4744 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4745 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4746 /* Record COO fields */ 4747 seqaij->coo_n = coo_n; 4748 seqaij->Atot = coo_n-nneg; /* Annz is seqaij->nz, so no need to record that again */ 4749 seqaij->jmap = jmap; /* of length nnz+1 */ 4750 seqaij->perm = perm; 4751 PetscFunctionReturn(0); 4752 } 4753 4754 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A,const PetscScalar v[],InsertMode imode) 4755 { 4756 Mat_SeqAIJ *aseq = (Mat_SeqAIJ*)A->data; 4757 PetscCount i,j,Annz = aseq->nz; 4758 PetscCount *perm = aseq->perm,*jmap = aseq->jmap; 4759 PetscScalar *Aa; 4760 4761 PetscFunctionBegin; 4762 PetscCall(MatSeqAIJGetArray(A,&Aa)); 4763 for (i=0; i<Annz; i++) { 4764 PetscScalar sum = 0.0; 4765 for (j=jmap[i]; j<jmap[i+1]; j++) sum += v[perm[j]]; 4766 Aa[i] = (imode == INSERT_VALUES? 0.0 : Aa[i]) + sum; 4767 } 4768 PetscCall(MatSeqAIJRestoreArray(A,&Aa)); 4769 PetscFunctionReturn(0); 4770 } 4771 4772 #if defined(PETSC_HAVE_CUDA) 4773 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat,MatType,MatReuse,Mat*); 4774 #endif 4775 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4776 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat,MatType,MatReuse,Mat*); 4777 #endif 4778 4779 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4780 { 4781 Mat_SeqAIJ *b; 4782 PetscMPIInt size; 4783 4784 PetscFunctionBegin; 4785 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 4786 PetscCheck(size <= 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4787 4788 PetscCall(PetscNewLog(B,&b)); 4789 4790 B->data = (void*)b; 4791 4792 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 4793 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4794 4795 b->row = NULL; 4796 b->col = NULL; 4797 b->icol = NULL; 4798 b->reallocs = 0; 4799 b->ignorezeroentries = PETSC_FALSE; 4800 b->roworiented = PETSC_TRUE; 4801 b->nonew = 0; 4802 b->diag = NULL; 4803 b->solve_work = NULL; 4804 B->spptr = NULL; 4805 b->saved_values = NULL; 4806 b->idiag = NULL; 4807 b->mdiag = NULL; 4808 b->ssor_work = NULL; 4809 b->omega = 1.0; 4810 b->fshift = 0.0; 4811 b->idiagvalid = PETSC_FALSE; 4812 b->ibdiagvalid = PETSC_FALSE; 4813 b->keepnonzeropattern = PETSC_FALSE; 4814 4815 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4816 4817 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4818 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ)); 4819 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ)); 4820 #endif 4821 4822 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ)); 4823 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ)); 4824 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ)); 4825 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ)); 4826 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ)); 4827 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM)); 4828 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL)); 4829 #if defined(PETSC_HAVE_MKL_SPARSE) 4830 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL)); 4831 #endif 4832 #if defined(PETSC_HAVE_CUDA) 4833 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4834 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4835 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",MatProductSetFromOptions_SeqAIJ)); 4836 #endif 4837 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4838 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijkokkos_C",MatConvert_SeqAIJ_SeqAIJKokkos)); 4839 #endif 4840 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL)); 4841 #if defined(PETSC_HAVE_ELEMENTAL) 4842 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental)); 4843 #endif 4844 #if defined(PETSC_HAVE_SCALAPACK) 4845 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK)); 4846 #endif 4847 #if defined(PETSC_HAVE_HYPRE) 4848 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE)); 4849 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4850 #endif 4851 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense)); 4852 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL)); 4853 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS)); 4854 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ)); 4855 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ)); 4856 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ)); 4857 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ)); 4858 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4859 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ)); 4860 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ)); 4861 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ)); 4862 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4863 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJKron_C",MatSeqAIJKron_SeqAIJ)); 4864 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_SeqAIJ)); 4865 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_SeqAIJ)); 4866 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4867 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4868 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4869 PetscFunctionReturn(0); 4870 } 4871 4872 /* 4873 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4874 */ 4875 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4876 { 4877 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data,*a = (Mat_SeqAIJ*)A->data; 4878 PetscInt m = A->rmap->n,i; 4879 4880 PetscFunctionBegin; 4881 PetscCheck(A->assembled || cpvalues == MAT_DO_NOT_COPY_VALUES,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot duplicate unassembled matrix"); 4882 4883 C->factortype = A->factortype; 4884 c->row = NULL; 4885 c->col = NULL; 4886 c->icol = NULL; 4887 c->reallocs = 0; 4888 4889 C->assembled = A->assembled; 4890 C->preallocated = A->preallocated; 4891 4892 if (A->preallocated) { 4893 PetscCall(PetscLayoutReference(A->rmap,&C->rmap)); 4894 PetscCall(PetscLayoutReference(A->cmap,&C->cmap)); 4895 4896 PetscCall(PetscMalloc1(m,&c->imax)); 4897 PetscCall(PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt))); 4898 PetscCall(PetscMalloc1(m,&c->ilen)); 4899 PetscCall(PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt))); 4900 PetscCall(PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt))); 4901 4902 /* allocate the matrix space */ 4903 if (mallocmatspace) { 4904 PetscCall(PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i)); 4905 PetscCall(PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt))); 4906 4907 c->singlemalloc = PETSC_TRUE; 4908 4909 PetscCall(PetscArraycpy(c->i,a->i,m+1)); 4910 if (m > 0) { 4911 PetscCall(PetscArraycpy(c->j,a->j,a->i[m])); 4912 if (cpvalues == MAT_COPY_VALUES) { 4913 const PetscScalar *aa; 4914 4915 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4916 PetscCall(PetscArraycpy(c->a,aa,a->i[m])); 4917 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4918 } else { 4919 PetscCall(PetscArrayzero(c->a,a->i[m])); 4920 } 4921 } 4922 } 4923 4924 c->ignorezeroentries = a->ignorezeroentries; 4925 c->roworiented = a->roworiented; 4926 c->nonew = a->nonew; 4927 if (a->diag) { 4928 PetscCall(PetscMalloc1(m+1,&c->diag)); 4929 PetscCall(PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt))); 4930 PetscCall(PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt))); 4931 } else c->diag = NULL; 4932 4933 c->solve_work = NULL; 4934 c->saved_values = NULL; 4935 c->idiag = NULL; 4936 c->ssor_work = NULL; 4937 c->keepnonzeropattern = a->keepnonzeropattern; 4938 c->free_a = PETSC_TRUE; 4939 c->free_ij = PETSC_TRUE; 4940 4941 c->rmax = a->rmax; 4942 c->nz = a->nz; 4943 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4944 4945 c->compressedrow.use = a->compressedrow.use; 4946 c->compressedrow.nrows = a->compressedrow.nrows; 4947 if (a->compressedrow.use) { 4948 i = a->compressedrow.nrows; 4949 PetscCall(PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex)); 4950 PetscCall(PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1)); 4951 PetscCall(PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i)); 4952 } else { 4953 c->compressedrow.use = PETSC_FALSE; 4954 c->compressedrow.i = NULL; 4955 c->compressedrow.rindex = NULL; 4956 } 4957 c->nonzerorowcnt = a->nonzerorowcnt; 4958 C->nonzerostate = A->nonzerostate; 4959 4960 PetscCall(MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C)); 4961 } 4962 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist)); 4963 PetscFunctionReturn(0); 4964 } 4965 4966 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4967 { 4968 PetscFunctionBegin; 4969 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 4970 PetscCall(MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 4971 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4972 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 4973 } 4974 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 4975 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE)); 4976 PetscFunctionReturn(0); 4977 } 4978 4979 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4980 { 4981 PetscBool isbinary, ishdf5; 4982 4983 PetscFunctionBegin; 4984 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4985 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4986 /* force binary viewer to load .info file if it has not yet done so */ 4987 PetscCall(PetscViewerSetUp(viewer)); 4988 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 4989 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5)); 4990 if (isbinary) { 4991 PetscCall(MatLoad_SeqAIJ_Binary(newMat,viewer)); 4992 } else if (ishdf5) { 4993 #if defined(PETSC_HAVE_HDF5) 4994 PetscCall(MatLoad_AIJ_HDF5(newMat,viewer)); 4995 #else 4996 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4997 #endif 4998 } else { 4999 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 5000 } 5001 PetscFunctionReturn(0); 5002 } 5003 5004 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 5005 { 5006 Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data; 5007 PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i; 5008 5009 PetscFunctionBegin; 5010 PetscCall(PetscViewerSetUp(viewer)); 5011 5012 /* read in matrix header */ 5013 PetscCall(PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT)); 5014 PetscCheck(header[0] == MAT_FILE_CLASSID,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 5015 M = header[1]; N = header[2]; nz = header[3]; 5016 PetscCheck(M >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%" PetscInt_FMT ") in file is negative",M); 5017 PetscCheck(N >= 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%" PetscInt_FMT ") in file is negative",N); 5018 PetscCheck(nz >= 0,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 5019 5020 /* set block sizes from the viewer's .info file */ 5021 PetscCall(MatLoad_Binary_BlockSizes(mat,viewer)); 5022 /* set local and global sizes if not set already */ 5023 if (mat->rmap->n < 0) mat->rmap->n = M; 5024 if (mat->cmap->n < 0) mat->cmap->n = N; 5025 if (mat->rmap->N < 0) mat->rmap->N = M; 5026 if (mat->cmap->N < 0) mat->cmap->N = N; 5027 PetscCall(PetscLayoutSetUp(mat->rmap)); 5028 PetscCall(PetscLayoutSetUp(mat->cmap)); 5029 5030 /* check if the matrix sizes are correct */ 5031 PetscCall(MatGetSize(mat,&rows,&cols)); 5032 PetscCheck(M == rows && N == cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5033 5034 /* read in row lengths */ 5035 PetscCall(PetscMalloc1(M,&rowlens)); 5036 PetscCall(PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT)); 5037 /* check if sum(rowlens) is same as nz */ 5038 sum = 0; for (i=0; i<M; i++) sum += rowlens[i]; 5039 PetscCheck(sum == nz,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT,nz,sum); 5040 /* preallocate and check sizes */ 5041 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens)); 5042 PetscCall(MatGetSize(mat,&rows,&cols)); 5043 PetscCheck(M == rows && N == cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5044 /* store row lengths */ 5045 PetscCall(PetscArraycpy(a->ilen,rowlens,M)); 5046 PetscCall(PetscFree(rowlens)); 5047 5048 /* fill in "i" row pointers */ 5049 a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i]; 5050 /* read in "j" column indices */ 5051 PetscCall(PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT)); 5052 /* read in "a" nonzero values */ 5053 PetscCall(PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR)); 5054 5055 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 5056 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 5057 PetscFunctionReturn(0); 5058 } 5059 5060 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 5061 { 5062 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 5063 const PetscScalar *aa,*ba; 5064 #if defined(PETSC_USE_COMPLEX) 5065 PetscInt k; 5066 #endif 5067 5068 PetscFunctionBegin; 5069 /* If the matrix dimensions are not equal,or no of nonzeros */ 5070 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 5071 *flg = PETSC_FALSE; 5072 PetscFunctionReturn(0); 5073 } 5074 5075 /* if the a->i are the same */ 5076 PetscCall(PetscArraycmp(a->i,b->i,A->rmap->n+1,flg)); 5077 if (!*flg) PetscFunctionReturn(0); 5078 5079 /* if a->j are the same */ 5080 PetscCall(PetscArraycmp(a->j,b->j,a->nz,flg)); 5081 if (!*flg) PetscFunctionReturn(0); 5082 5083 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 5084 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 5085 /* if a->a are the same */ 5086 #if defined(PETSC_USE_COMPLEX) 5087 for (k=0; k<a->nz; k++) { 5088 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5089 *flg = PETSC_FALSE; 5090 PetscFunctionReturn(0); 5091 } 5092 } 5093 #else 5094 PetscCall(PetscArraycmp(aa,ba,a->nz,flg)); 5095 #endif 5096 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 5097 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 5098 PetscFunctionReturn(0); 5099 } 5100 5101 /*@ 5102 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 5103 provided by the user. 5104 5105 Collective 5106 5107 Input Parameters: 5108 + comm - must be an MPI communicator of size 1 5109 . m - number of rows 5110 . n - number of columns 5111 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5112 . j - column indices 5113 - a - matrix values 5114 5115 Output Parameter: 5116 . mat - the matrix 5117 5118 Level: intermediate 5119 5120 Notes: 5121 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5122 once the matrix is destroyed and not before 5123 5124 You cannot set new nonzero locations into this matrix, that will generate an error. 5125 5126 The i and j indices are 0 based 5127 5128 The format which is used for the sparse matrix input, is equivalent to a 5129 row-major ordering.. i.e for the following matrix, the input data expected is 5130 as shown 5131 5132 $ 1 0 0 5133 $ 2 0 3 5134 $ 4 5 6 5135 $ 5136 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5137 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5138 $ v = {1,2,3,4,5,6} [size = 6] 5139 5140 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateMPIAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()` 5141 5142 @*/ 5143 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 5144 { 5145 PetscInt ii; 5146 Mat_SeqAIJ *aij; 5147 PetscInt jj; 5148 5149 PetscFunctionBegin; 5150 PetscCheck(m <= 0 || i[0] == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 5151 PetscCall(MatCreate(comm,mat)); 5152 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5153 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5154 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5155 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL)); 5156 aij = (Mat_SeqAIJ*)(*mat)->data; 5157 PetscCall(PetscMalloc1(m,&aij->imax)); 5158 PetscCall(PetscMalloc1(m,&aij->ilen)); 5159 5160 aij->i = i; 5161 aij->j = j; 5162 aij->a = a; 5163 aij->singlemalloc = PETSC_FALSE; 5164 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5165 aij->free_a = PETSC_FALSE; 5166 aij->free_ij = PETSC_FALSE; 5167 5168 for (ii=0,aij->nonzerorowcnt=0,aij->rmax=0; ii<m; ii++) { 5169 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 5170 if (PetscDefined(USE_DEBUG)) { 5171 PetscCheck(i[ii+1] - i[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT,ii,i[ii+1] - i[ii]); 5172 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 5173 PetscCheck(j[jj] >= j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is not sorted",jj-i[ii],j[jj],ii); 5174 PetscCheck(j[jj] != j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is identical to previous entry",jj-i[ii],j[jj],ii); 5175 } 5176 } 5177 } 5178 if (PetscDefined(USE_DEBUG)) { 5179 for (ii=0; ii<aij->i[m]; ii++) { 5180 PetscCheck(j[ii] >= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5181 PetscCheck(j[ii] <= n - 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5182 } 5183 } 5184 5185 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5186 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5187 PetscFunctionReturn(0); 5188 } 5189 5190 /*@ 5191 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 5192 provided by the user. 5193 5194 Collective 5195 5196 Input Parameters: 5197 + comm - must be an MPI communicator of size 1 5198 . m - number of rows 5199 . n - number of columns 5200 . i - row indices 5201 . j - column indices 5202 . a - matrix values 5203 . nz - number of nonzeros 5204 - idx - if the i and j indices start with 1 use PETSC_TRUE otherwise use PETSC_FALSE 5205 5206 Output Parameter: 5207 . mat - the matrix 5208 5209 Level: intermediate 5210 5211 Example: 5212 For the following matrix, the input data expected is as shown (using 0 based indexing) 5213 .vb 5214 1 0 0 5215 2 0 3 5216 4 5 6 5217 5218 i = {0,1,1,2,2,2} 5219 j = {0,0,2,0,1,2} 5220 v = {1,2,3,4,5,6} 5221 .ve 5222 5223 .seealso: `MatCreate()`, `MatCreateAIJ()`, `MatCreateSeqAIJ()`, `MatCreateSeqAIJWithArrays()`, `MatMPIAIJSetPreallocationCSR()`, `MatSetValuesCOO()` 5224 5225 @*/ 5226 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 5227 { 5228 PetscInt ii, *nnz, one = 1,row,col; 5229 5230 PetscFunctionBegin; 5231 PetscCall(PetscCalloc1(m,&nnz)); 5232 for (ii = 0; ii < nz; ii++) { 5233 nnz[i[ii] - !!idx] += 1; 5234 } 5235 PetscCall(MatCreate(comm,mat)); 5236 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5237 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5238 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz)); 5239 for (ii = 0; ii < nz; ii++) { 5240 if (idx) { 5241 row = i[ii] - 1; 5242 col = j[ii] - 1; 5243 } else { 5244 row = i[ii]; 5245 col = j[ii]; 5246 } 5247 PetscCall(MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES)); 5248 } 5249 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5250 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5251 PetscCall(PetscFree(nnz)); 5252 PetscFunctionReturn(0); 5253 } 5254 5255 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5256 { 5257 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 5258 5259 PetscFunctionBegin; 5260 a->idiagvalid = PETSC_FALSE; 5261 a->ibdiagvalid = PETSC_FALSE; 5262 5263 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5264 PetscFunctionReturn(0); 5265 } 5266 5267 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 5268 { 5269 PetscMPIInt size; 5270 5271 PetscFunctionBegin; 5272 PetscCallMPI(MPI_Comm_size(comm,&size)); 5273 if (size == 1) { 5274 if (scall == MAT_INITIAL_MATRIX) { 5275 PetscCall(MatDuplicate(inmat,MAT_COPY_VALUES,outmat)); 5276 } else { 5277 PetscCall(MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN)); 5278 } 5279 } else { 5280 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat)); 5281 } 5282 PetscFunctionReturn(0); 5283 } 5284 5285 /* 5286 Permute A into C's *local* index space using rowemb,colemb. 5287 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5288 of [0,m), colemb is in [0,n). 5289 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5290 */ 5291 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 5292 { 5293 /* If making this function public, change the error returned in this function away from _PLIB. */ 5294 Mat_SeqAIJ *Baij; 5295 PetscBool seqaij; 5296 PetscInt m,n,*nz,i,j,count; 5297 PetscScalar v; 5298 const PetscInt *rowindices,*colindices; 5299 5300 PetscFunctionBegin; 5301 if (!B) PetscFunctionReturn(0); 5302 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5303 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij)); 5304 PetscCheck(seqaij,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 5305 if (rowemb) { 5306 PetscCall(ISGetLocalSize(rowemb,&m)); 5307 PetscCheck(m == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Row IS of size %" PetscInt_FMT " is incompatible with matrix row size %" PetscInt_FMT,m,B->rmap->n); 5308 } else { 5309 PetscCheck(C->rmap->n == B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 5310 } 5311 if (colemb) { 5312 PetscCall(ISGetLocalSize(colemb,&n)); 5313 PetscCheck(n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diag col IS of size %" PetscInt_FMT " is incompatible with input matrix col size %" PetscInt_FMT,n,B->cmap->n); 5314 } else { 5315 PetscCheck(C->cmap->n == B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 5316 } 5317 5318 Baij = (Mat_SeqAIJ*)(B->data); 5319 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5320 PetscCall(PetscMalloc1(B->rmap->n,&nz)); 5321 for (i=0; i<B->rmap->n; i++) { 5322 nz[i] = Baij->i[i+1] - Baij->i[i]; 5323 } 5324 PetscCall(MatSeqAIJSetPreallocation(C,0,nz)); 5325 PetscCall(PetscFree(nz)); 5326 } 5327 if (pattern == SUBSET_NONZERO_PATTERN) { 5328 PetscCall(MatZeroEntries(C)); 5329 } 5330 count = 0; 5331 rowindices = NULL; 5332 colindices = NULL; 5333 if (rowemb) { 5334 PetscCall(ISGetIndices(rowemb,&rowindices)); 5335 } 5336 if (colemb) { 5337 PetscCall(ISGetIndices(colemb,&colindices)); 5338 } 5339 for (i=0; i<B->rmap->n; i++) { 5340 PetscInt row; 5341 row = i; 5342 if (rowindices) row = rowindices[i]; 5343 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 5344 PetscInt col; 5345 col = Baij->j[count]; 5346 if (colindices) col = colindices[col]; 5347 v = Baij->a[count]; 5348 PetscCall(MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES)); 5349 ++count; 5350 } 5351 } 5352 /* FIXME: set C's nonzerostate correctly. */ 5353 /* Assembly for C is necessary. */ 5354 C->preallocated = PETSC_TRUE; 5355 C->assembled = PETSC_TRUE; 5356 C->was_assembled = PETSC_FALSE; 5357 PetscFunctionReturn(0); 5358 } 5359 5360 PetscFunctionList MatSeqAIJList = NULL; 5361 5362 /*@C 5363 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 5364 5365 Collective on Mat 5366 5367 Input Parameters: 5368 + mat - the matrix object 5369 - matype - matrix type 5370 5371 Options Database Key: 5372 . -mat_seqai_type <method> - for example seqaijcrl 5373 5374 Level: intermediate 5375 5376 .seealso: `PCSetType()`, `VecSetType()`, `MatCreate()`, `MatType`, `Mat` 5377 @*/ 5378 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5379 { 5380 PetscBool sametype; 5381 PetscErrorCode (*r)(Mat,MatType,MatReuse,Mat*); 5382 5383 PetscFunctionBegin; 5384 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 5385 PetscCall(PetscObjectTypeCompare((PetscObject)mat,matype,&sametype)); 5386 if (sametype) PetscFunctionReturn(0); 5387 5388 PetscCall(PetscFunctionListFind(MatSeqAIJList,matype,&r)); 5389 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 5390 PetscCall((*r)(mat,matype,MAT_INPLACE_MATRIX,&mat)); 5391 PetscFunctionReturn(0); 5392 } 5393 5394 /*@C 5395 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 5396 5397 Not Collective 5398 5399 Input Parameters: 5400 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 5401 - function - routine to convert to subtype 5402 5403 Notes: 5404 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 5405 5406 Then, your matrix can be chosen with the procedural interface at runtime via the option 5407 $ -mat_seqaij_type my_mat 5408 5409 Level: advanced 5410 5411 .seealso: `MatSeqAIJRegisterAll()` 5412 5413 Level: advanced 5414 @*/ 5415 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 5416 { 5417 PetscFunctionBegin; 5418 PetscCall(MatInitializePackage()); 5419 PetscCall(PetscFunctionListAdd(&MatSeqAIJList,sname,function)); 5420 PetscFunctionReturn(0); 5421 } 5422 5423 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5424 5425 /*@C 5426 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 5427 5428 Not Collective 5429 5430 Level: advanced 5431 5432 .seealso: `MatRegisterAll()`, `MatSeqAIJRegister()` 5433 @*/ 5434 PetscErrorCode MatSeqAIJRegisterAll(void) 5435 { 5436 PetscFunctionBegin; 5437 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5438 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5439 5440 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5441 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5442 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5443 #if defined(PETSC_HAVE_MKL_SPARSE) 5444 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5445 #endif 5446 #if defined(PETSC_HAVE_CUDA) 5447 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5448 #endif 5449 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5450 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5451 #endif 5452 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5453 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5454 #endif 5455 PetscFunctionReturn(0); 5456 } 5457 5458 /* 5459 Special version for direct calls from Fortran 5460 */ 5461 #include <petsc/private/fortranimpl.h> 5462 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5463 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5464 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5465 #define matsetvaluesseqaij_ matsetvaluesseqaij 5466 #endif 5467 5468 /* Change these macros so can be used in void function */ 5469 5470 /* Change these macros so can be used in void function */ 5471 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5472 #undef PetscCall 5473 #define PetscCall(...) do { \ 5474 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5475 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5476 *_ierr = PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr_msv_mpiaij,PETSC_ERROR_REPEAT," "); \ 5477 return; \ 5478 } \ 5479 } while (0) 5480 5481 #undef SETERRQ 5482 #define SETERRQ(comm,ierr,...) do { \ 5483 *_ierr = PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__); \ 5484 return; \ 5485 } while (0) 5486 5487 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 5488 { 5489 Mat A = *AA; 5490 PetscInt m = *mm, n = *nn; 5491 InsertMode is = *isis; 5492 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 5493 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 5494 PetscInt *imax,*ai,*ailen; 5495 PetscInt *aj,nonew = a->nonew,lastcol = -1; 5496 MatScalar *ap,value,*aa; 5497 PetscBool ignorezeroentries = a->ignorezeroentries; 5498 PetscBool roworiented = a->roworiented; 5499 5500 PetscFunctionBegin; 5501 MatCheckPreallocated(A,1); 5502 imax = a->imax; 5503 ai = a->i; 5504 ailen = a->ilen; 5505 aj = a->j; 5506 aa = a->a; 5507 5508 for (k=0; k<m; k++) { /* loop over added rows */ 5509 row = im[k]; 5510 if (row < 0) continue; 5511 PetscCheck(row < A->rmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 5512 rp = aj + ai[row]; ap = aa + ai[row]; 5513 rmax = imax[row]; nrow = ailen[row]; 5514 low = 0; 5515 high = nrow; 5516 for (l=0; l<n; l++) { /* loop over added columns */ 5517 if (in[l] < 0) continue; 5518 PetscCheck(in[l] < A->cmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 5519 col = in[l]; 5520 if (roworiented) value = v[l + k*n]; 5521 else value = v[k + l*m]; 5522 5523 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5524 5525 if (col <= lastcol) low = 0; 5526 else high = nrow; 5527 lastcol = col; 5528 while (high-low > 5) { 5529 t = (low+high)/2; 5530 if (rp[t] > col) high = t; 5531 else low = t; 5532 } 5533 for (i=low; i<high; i++) { 5534 if (rp[i] > col) break; 5535 if (rp[i] == col) { 5536 if (is == ADD_VALUES) ap[i] += value; 5537 else ap[i] = value; 5538 goto noinsert; 5539 } 5540 } 5541 if (value == 0.0 && ignorezeroentries) goto noinsert; 5542 if (nonew == 1) goto noinsert; 5543 PetscCheck(nonew != -1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 5544 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 5545 N = nrow++ - 1; a->nz++; high++; 5546 /* shift up all the later entries in this row */ 5547 for (ii=N; ii>=i; ii--) { 5548 rp[ii+1] = rp[ii]; 5549 ap[ii+1] = ap[ii]; 5550 } 5551 rp[i] = col; 5552 ap[i] = value; 5553 A->nonzerostate++; 5554 noinsert:; 5555 low = i + 1; 5556 } 5557 ailen[row] = nrow; 5558 } 5559 PetscFunctionReturnVoid(); 5560 } 5561 /* Undefining these here since they were redefined from their original definition above! No 5562 * other PETSc functions should be defined past this point, as it is impossible to recover the 5563 * original definitions */ 5564 #undef PetscCall 5565 #undef SETERRQ 5566