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