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