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