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