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