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