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