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 = NULL; 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,NULL);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=NULL,*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 = NULL; 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,NULL);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 NULL, 3383 NULL, 3384 NULL, 3385 /* 10*/ NULL, 3386 MatLUFactor_SeqAIJ, 3387 NULL, 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*/ NULL, 3396 MatAssemblyEnd_SeqAIJ, 3397 MatSetOption_SeqAIJ, 3398 MatZeroEntries_SeqAIJ, 3399 /* 24*/ MatZeroRows_SeqAIJ, 3400 NULL, 3401 NULL, 3402 NULL, 3403 NULL, 3404 /* 29*/ MatSetUp_SeqAIJ, 3405 NULL, 3406 NULL, 3407 NULL, 3408 NULL, 3409 /* 34*/ MatDuplicate_SeqAIJ, 3410 NULL, 3411 NULL, 3412 MatILUFactor_SeqAIJ, 3413 NULL, 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 NULL, 3431 NULL, 3432 MatPermute_SeqAIJ, 3433 NULL, 3434 /* 59*/ NULL, 3435 MatDestroy_SeqAIJ, 3436 MatView_SeqAIJ, 3437 NULL, 3438 NULL, 3439 /* 64*/ NULL, 3440 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3441 NULL, 3442 NULL, 3443 NULL, 3444 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3445 MatGetRowMinAbs_SeqAIJ, 3446 NULL, 3447 NULL, 3448 NULL, 3449 /* 74*/ NULL, 3450 MatFDColoringApply_AIJ, 3451 NULL, 3452 NULL, 3453 NULL, 3454 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3455 NULL, 3456 NULL, 3457 NULL, 3458 MatLoad_SeqAIJ, 3459 /* 84*/ MatIsSymmetric_SeqAIJ, 3460 MatIsHermitian_SeqAIJ, 3461 NULL, 3462 NULL, 3463 NULL, 3464 /* 89*/ NULL, 3465 NULL, 3466 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3467 NULL, 3468 NULL, 3469 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3470 NULL, 3471 NULL, 3472 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3473 NULL, 3474 /* 99*/ MatProductSetFromOptions_SeqAIJ, 3475 NULL, 3476 NULL, 3477 MatConjugate_SeqAIJ, 3478 NULL, 3479 /*104*/ MatSetValuesRow_SeqAIJ, 3480 MatRealPart_SeqAIJ, 3481 MatImaginaryPart_SeqAIJ, 3482 NULL, 3483 NULL, 3484 /*109*/ MatMatSolve_SeqAIJ, 3485 NULL, 3486 MatGetRowMin_SeqAIJ, 3487 NULL, 3488 MatMissingDiagonal_SeqAIJ, 3489 /*114*/ NULL, 3490 NULL, 3491 NULL, 3492 NULL, 3493 NULL, 3494 /*119*/ NULL, 3495 NULL, 3496 NULL, 3497 NULL, 3498 MatGetMultiProcBlock_SeqAIJ, 3499 /*124*/ MatFindNonzeroRows_SeqAIJ, 3500 MatGetColumnNorms_SeqAIJ, 3501 MatInvertBlockDiagonal_SeqAIJ, 3502 MatInvertVariableBlockDiagonal_SeqAIJ, 3503 NULL, 3504 /*129*/ NULL, 3505 NULL, 3506 NULL, 3507 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3508 MatTransposeColoringCreate_SeqAIJ, 3509 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3510 MatTransColoringApplyDenToSp_SeqAIJ, 3511 NULL, 3512 NULL, 3513 MatRARtNumeric_SeqAIJ_SeqAIJ, 3514 /*139*/NULL, 3515 NULL, 3516 NULL, 3517 MatFDColoringSetUp_SeqXAIJ, 3518 MatFindOffBlockDiagonalEntries_SeqAIJ, 3519 MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3520 /*145*/MatDestroySubMatrices_SeqAIJ, 3521 NULL, 3522 NULL 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 Notes: 4034 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 4035 4036 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4037 structure will be the union of all the previous nonzero structures. 4038 4039 Developer Notes: 4040 An optimization could be added to the implementation where it checks if the i, and j are identical to the current i and j and 4041 then just copies the v values directly with PetscMemcpy(). 4042 4043 This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them. 4044 4045 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ, MatResetPreallocation() 4046 @*/ 4047 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 4048 { 4049 PetscErrorCode ierr; 4050 4051 PetscFunctionBegin; 4052 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 4053 PetscValidType(B,1); 4054 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 4055 PetscFunctionReturn(0); 4056 } 4057 4058 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 4059 { 4060 PetscInt i; 4061 PetscInt m,n; 4062 PetscInt nz; 4063 PetscInt *nnz; 4064 PetscErrorCode ierr; 4065 4066 PetscFunctionBegin; 4067 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 4068 4069 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 4070 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 4071 4072 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 4073 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 4074 for (i = 0; i < m; i++) { 4075 nz = Ii[i+1]- Ii[i]; 4076 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 4077 nnz[i] = nz; 4078 } 4079 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 4080 ierr = PetscFree(nnz);CHKERRQ(ierr); 4081 4082 for (i = 0; i < m; i++) { 4083 ierr = MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES);CHKERRQ(ierr); 4084 } 4085 4086 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4087 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4088 4089 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 4090 PetscFunctionReturn(0); 4091 } 4092 4093 #include <../src/mat/impls/dense/seq/dense.h> 4094 #include <petsc/private/kernels/petscaxpy.h> 4095 4096 /* 4097 Computes (B'*A')' since computing B*A directly is untenable 4098 4099 n p p 4100 [ ] [ ] [ ] 4101 m [ A ] * n [ B ] = m [ C ] 4102 [ ] [ ] [ ] 4103 4104 */ 4105 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 4106 { 4107 PetscErrorCode ierr; 4108 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 4109 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 4110 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 4111 PetscInt i,j,n,m,q,p; 4112 const PetscInt *ii,*idx; 4113 const PetscScalar *b,*a,*a_q; 4114 PetscScalar *c,*c_q; 4115 PetscInt clda = sub_c->lda; 4116 PetscInt alda = sub_a->lda; 4117 4118 PetscFunctionBegin; 4119 m = A->rmap->n; 4120 n = A->cmap->n; 4121 p = B->cmap->n; 4122 a = sub_a->v; 4123 b = sub_b->a; 4124 c = sub_c->v; 4125 if (clda == m) { 4126 ierr = PetscArrayzero(c,m*p);CHKERRQ(ierr); 4127 } else { 4128 for (j=0;j<p;j++) 4129 for (i=0;i<m;i++) 4130 c[j*clda + i] = 0.0; 4131 } 4132 ii = sub_b->i; 4133 idx = sub_b->j; 4134 for (i=0; i<n; i++) { 4135 q = ii[i+1] - ii[i]; 4136 while (q-->0) { 4137 c_q = c + clda*(*idx); 4138 a_q = a + alda*i; 4139 PetscKernelAXPY(c_q,*b,a_q,m); 4140 idx++; 4141 b++; 4142 } 4143 } 4144 PetscFunctionReturn(0); 4145 } 4146 4147 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C) 4148 { 4149 PetscErrorCode ierr; 4150 PetscInt m=A->rmap->n,n=B->cmap->n; 4151 PetscBool cisdense; 4152 4153 PetscFunctionBegin; 4154 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); 4155 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 4156 ierr = MatSetBlockSizesFromMats(C,A,B);CHKERRQ(ierr); 4157 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 4158 if (!cisdense) { 4159 ierr = MatSetType(C,MATDENSE);CHKERRQ(ierr); 4160 } 4161 ierr = MatSetUp(C);CHKERRQ(ierr); 4162 4163 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4164 PetscFunctionReturn(0); 4165 } 4166 4167 /* ----------------------------------------------------------------*/ 4168 /*MC 4169 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4170 based on compressed sparse row format. 4171 4172 Options Database Keys: 4173 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4174 4175 Level: beginner 4176 4177 Notes: 4178 MatSetValues() may be called for this matrix type with a NULL argument for the numerical values, 4179 in this case the values associated with the rows and columns one passes in are set to zero 4180 in the matrix 4181 4182 MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no 4183 space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored 4184 4185 Developer Notes: 4186 It would be nice if all matrix formats supported passing NULL in for the numerical values 4187 4188 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 4189 M*/ 4190 4191 /*MC 4192 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4193 4194 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4195 and MATMPIAIJ otherwise. As a result, for single process communicators, 4196 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation() is supported 4197 for communicators controlling multiple processes. It is recommended that you call both of 4198 the above preallocation routines for simplicity. 4199 4200 Options Database Keys: 4201 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4202 4203 Developer Notes: 4204 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4205 enough exist. 4206 4207 Level: beginner 4208 4209 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 4210 M*/ 4211 4212 /*MC 4213 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4214 4215 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4216 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4217 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4218 for communicators controlling multiple processes. It is recommended that you call both of 4219 the above preallocation routines for simplicity. 4220 4221 Options Database Keys: 4222 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4223 4224 Level: beginner 4225 4226 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 4227 M*/ 4228 4229 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4230 #if defined(PETSC_HAVE_ELEMENTAL) 4231 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4232 #endif 4233 #if defined(PETSC_HAVE_SCALAPACK) 4234 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*); 4235 #endif 4236 #if defined(PETSC_HAVE_HYPRE) 4237 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4238 #endif 4239 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 4240 4241 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4242 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4243 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4244 4245 /*@C 4246 MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored 4247 4248 Not Collective 4249 4250 Input Parameter: 4251 . mat - a MATSEQAIJ matrix 4252 4253 Output Parameter: 4254 . array - pointer to the data 4255 4256 Level: intermediate 4257 4258 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4259 @*/ 4260 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4261 { 4262 PetscErrorCode ierr; 4263 4264 PetscFunctionBegin; 4265 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4266 PetscFunctionReturn(0); 4267 } 4268 4269 /*@C 4270 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored 4271 4272 Not Collective 4273 4274 Input Parameter: 4275 . mat - a MATSEQAIJ matrix 4276 4277 Output Parameter: 4278 . array - pointer to the data 4279 4280 Level: intermediate 4281 4282 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead() 4283 @*/ 4284 PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array) 4285 { 4286 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL) 4287 PetscOffloadMask oval; 4288 #endif 4289 PetscErrorCode ierr; 4290 4291 PetscFunctionBegin; 4292 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL) 4293 oval = A->offloadmask; 4294 #endif 4295 ierr = MatSeqAIJGetArray(A,(PetscScalar**)array);CHKERRQ(ierr); 4296 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL) 4297 if (oval == PETSC_OFFLOAD_GPU || oval == PETSC_OFFLOAD_BOTH) A->offloadmask = PETSC_OFFLOAD_BOTH; 4298 #endif 4299 PetscFunctionReturn(0); 4300 } 4301 4302 /*@C 4303 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4304 4305 Not Collective 4306 4307 Input Parameter: 4308 . mat - a MATSEQAIJ matrix 4309 4310 Output Parameter: 4311 . array - pointer to the data 4312 4313 Level: intermediate 4314 4315 .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead() 4316 @*/ 4317 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array) 4318 { 4319 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL) 4320 PetscOffloadMask oval; 4321 #endif 4322 PetscErrorCode ierr; 4323 4324 PetscFunctionBegin; 4325 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL) 4326 oval = A->offloadmask; 4327 #endif 4328 ierr = MatSeqAIJRestoreArray(A,(PetscScalar**)array);CHKERRQ(ierr); 4329 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_VIENNACL) 4330 A->offloadmask = oval; 4331 #endif 4332 PetscFunctionReturn(0); 4333 } 4334 4335 /*@C 4336 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4337 4338 Not Collective 4339 4340 Input Parameter: 4341 . mat - a MATSEQAIJ matrix 4342 4343 Output Parameter: 4344 . nz - the maximum number of nonzeros in any row 4345 4346 Level: intermediate 4347 4348 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4349 @*/ 4350 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4351 { 4352 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4353 4354 PetscFunctionBegin; 4355 *nz = aij->rmax; 4356 PetscFunctionReturn(0); 4357 } 4358 4359 /*@C 4360 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4361 4362 Not Collective 4363 4364 Input Parameters: 4365 + mat - a MATSEQAIJ matrix 4366 - array - pointer to the data 4367 4368 Level: intermediate 4369 4370 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4371 @*/ 4372 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4373 { 4374 PetscErrorCode ierr; 4375 4376 PetscFunctionBegin; 4377 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4378 PetscFunctionReturn(0); 4379 } 4380 4381 #if defined(PETSC_HAVE_CUDA) 4382 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat); 4383 #endif 4384 4385 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4386 { 4387 Mat_SeqAIJ *b; 4388 PetscErrorCode ierr; 4389 PetscMPIInt size; 4390 4391 PetscFunctionBegin; 4392 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4393 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4394 4395 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4396 4397 B->data = (void*)b; 4398 4399 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4400 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4401 4402 b->row = NULL; 4403 b->col = NULL; 4404 b->icol = NULL; 4405 b->reallocs = 0; 4406 b->ignorezeroentries = PETSC_FALSE; 4407 b->roworiented = PETSC_TRUE; 4408 b->nonew = 0; 4409 b->diag = NULL; 4410 b->solve_work = NULL; 4411 B->spptr = NULL; 4412 b->saved_values = NULL; 4413 b->idiag = NULL; 4414 b->mdiag = NULL; 4415 b->ssor_work = NULL; 4416 b->omega = 1.0; 4417 b->fshift = 0.0; 4418 b->idiagvalid = PETSC_FALSE; 4419 b->ibdiagvalid = PETSC_FALSE; 4420 b->keepnonzeropattern = PETSC_FALSE; 4421 4422 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4423 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4424 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4425 4426 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4427 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4428 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4429 #endif 4430 4431 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4432 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4433 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4434 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4435 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4436 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4437 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL);CHKERRQ(ierr); 4438 #if defined(PETSC_HAVE_MKL_SPARSE) 4439 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4440 #endif 4441 #if defined(PETSC_HAVE_CUDA) 4442 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE);CHKERRQ(ierr); 4443 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ);CHKERRQ(ierr); 4444 #endif 4445 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4446 #if defined(PETSC_HAVE_ELEMENTAL) 4447 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4448 #endif 4449 #if defined(PETSC_HAVE_SCALAPACK) 4450 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK);CHKERRQ(ierr); 4451 #endif 4452 #if defined(PETSC_HAVE_HYPRE) 4453 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 4454 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ);CHKERRQ(ierr); 4455 #endif 4456 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 4457 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);CHKERRQ(ierr); 4458 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS);CHKERRQ(ierr); 4459 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4460 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4461 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4462 ierr = PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);CHKERRQ(ierr); 4463 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4464 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4465 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ);CHKERRQ(ierr); 4466 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ);CHKERRQ(ierr); 4467 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ);CHKERRQ(ierr); 4468 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4469 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4470 ierr = MatSeqAIJSetTypeFromOptions(B);CHKERRQ(ierr); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4471 PetscFunctionReturn(0); 4472 } 4473 4474 /* 4475 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4476 */ 4477 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4478 { 4479 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4480 PetscErrorCode ierr; 4481 PetscInt m = A->rmap->n,i; 4482 4483 PetscFunctionBegin; 4484 c = (Mat_SeqAIJ*)C->data; 4485 4486 C->factortype = A->factortype; 4487 c->row = NULL; 4488 c->col = NULL; 4489 c->icol = NULL; 4490 c->reallocs = 0; 4491 4492 C->assembled = PETSC_TRUE; 4493 4494 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4495 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4496 4497 ierr = PetscMalloc1(m,&c->imax);CHKERRQ(ierr); 4498 ierr = PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt));CHKERRQ(ierr); 4499 ierr = PetscMalloc1(m,&c->ilen);CHKERRQ(ierr); 4500 ierr = PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt));CHKERRQ(ierr); 4501 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4502 4503 /* allocate the matrix space */ 4504 if (mallocmatspace) { 4505 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4506 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4507 4508 c->singlemalloc = PETSC_TRUE; 4509 4510 ierr = PetscArraycpy(c->i,a->i,m+1);CHKERRQ(ierr); 4511 if (m > 0) { 4512 ierr = PetscArraycpy(c->j,a->j,a->i[m]);CHKERRQ(ierr); 4513 if (cpvalues == MAT_COPY_VALUES) { 4514 ierr = PetscArraycpy(c->a,a->a,a->i[m]);CHKERRQ(ierr); 4515 } else { 4516 ierr = PetscArrayzero(c->a,a->i[m]);CHKERRQ(ierr); 4517 } 4518 } 4519 } 4520 4521 c->ignorezeroentries = a->ignorezeroentries; 4522 c->roworiented = a->roworiented; 4523 c->nonew = a->nonew; 4524 if (a->diag) { 4525 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4526 ierr = PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt));CHKERRQ(ierr); 4527 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4528 } else c->diag = NULL; 4529 4530 c->solve_work = NULL; 4531 c->saved_values = NULL; 4532 c->idiag = NULL; 4533 c->ssor_work = NULL; 4534 c->keepnonzeropattern = a->keepnonzeropattern; 4535 c->free_a = PETSC_TRUE; 4536 c->free_ij = PETSC_TRUE; 4537 4538 c->rmax = a->rmax; 4539 c->nz = a->nz; 4540 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4541 C->preallocated = PETSC_TRUE; 4542 4543 c->compressedrow.use = a->compressedrow.use; 4544 c->compressedrow.nrows = a->compressedrow.nrows; 4545 if (a->compressedrow.use) { 4546 i = a->compressedrow.nrows; 4547 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4548 ierr = PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1);CHKERRQ(ierr); 4549 ierr = PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i);CHKERRQ(ierr); 4550 } else { 4551 c->compressedrow.use = PETSC_FALSE; 4552 c->compressedrow.i = NULL; 4553 c->compressedrow.rindex = NULL; 4554 } 4555 c->nonzerorowcnt = a->nonzerorowcnt; 4556 C->nonzerostate = A->nonzerostate; 4557 4558 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4559 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4560 PetscFunctionReturn(0); 4561 } 4562 4563 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4564 { 4565 PetscErrorCode ierr; 4566 4567 PetscFunctionBegin; 4568 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4569 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4570 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4571 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4572 } 4573 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4574 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4575 PetscFunctionReturn(0); 4576 } 4577 4578 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4579 { 4580 PetscBool isbinary, ishdf5; 4581 PetscErrorCode ierr; 4582 4583 PetscFunctionBegin; 4584 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4585 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4586 /* force binary viewer to load .info file if it has not yet done so */ 4587 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4588 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 4589 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 4590 if (isbinary) { 4591 ierr = MatLoad_SeqAIJ_Binary(newMat,viewer);CHKERRQ(ierr); 4592 } else if (ishdf5) { 4593 #if defined(PETSC_HAVE_HDF5) 4594 ierr = MatLoad_AIJ_HDF5(newMat,viewer);CHKERRQ(ierr); 4595 #else 4596 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4597 #endif 4598 } else { 4599 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); 4600 } 4601 PetscFunctionReturn(0); 4602 } 4603 4604 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 4605 { 4606 Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data; 4607 PetscErrorCode ierr; 4608 PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i; 4609 4610 PetscFunctionBegin; 4611 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4612 4613 /* read in matrix header */ 4614 ierr = PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);CHKERRQ(ierr); 4615 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 4616 M = header[1]; N = header[2]; nz = header[3]; 4617 if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M); 4618 if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N); 4619 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 4620 4621 /* set block sizes from the viewer's .info file */ 4622 ierr = MatLoad_Binary_BlockSizes(mat,viewer);CHKERRQ(ierr); 4623 /* set local and global sizes if not set already */ 4624 if (mat->rmap->n < 0) mat->rmap->n = M; 4625 if (mat->cmap->n < 0) mat->cmap->n = N; 4626 if (mat->rmap->N < 0) mat->rmap->N = M; 4627 if (mat->cmap->N < 0) mat->cmap->N = N; 4628 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 4629 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 4630 4631 /* check if the matrix sizes are correct */ 4632 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 4633 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); 4634 4635 /* read in row lengths */ 4636 ierr = PetscMalloc1(M,&rowlens);CHKERRQ(ierr); 4637 ierr = PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT);CHKERRQ(ierr); 4638 /* check if sum(rowlens) is same as nz */ 4639 sum = 0; for (i=0; i<M; i++) sum += rowlens[i]; 4640 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); 4641 /* preallocate and check sizes */ 4642 ierr = MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens);CHKERRQ(ierr); 4643 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 4644 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); 4645 /* store row lengths */ 4646 ierr = PetscArraycpy(a->ilen,rowlens,M);CHKERRQ(ierr); 4647 ierr = PetscFree(rowlens);CHKERRQ(ierr); 4648 4649 /* fill in "i" row pointers */ 4650 a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i]; 4651 /* read in "j" column indices */ 4652 ierr = PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT);CHKERRQ(ierr); 4653 /* read in "a" nonzero values */ 4654 ierr = PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR);CHKERRQ(ierr); 4655 4656 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4657 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4658 PetscFunctionReturn(0); 4659 } 4660 4661 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4662 { 4663 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4664 PetscErrorCode ierr; 4665 #if defined(PETSC_USE_COMPLEX) 4666 PetscInt k; 4667 #endif 4668 4669 PetscFunctionBegin; 4670 /* If the matrix dimensions are not equal,or no of nonzeros */ 4671 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4672 *flg = PETSC_FALSE; 4673 PetscFunctionReturn(0); 4674 } 4675 4676 /* if the a->i are the same */ 4677 ierr = PetscArraycmp(a->i,b->i,A->rmap->n+1,flg);CHKERRQ(ierr); 4678 if (!*flg) PetscFunctionReturn(0); 4679 4680 /* if a->j are the same */ 4681 ierr = PetscArraycmp(a->j,b->j,a->nz,flg);CHKERRQ(ierr); 4682 if (!*flg) PetscFunctionReturn(0); 4683 4684 /* if a->a are the same */ 4685 #if defined(PETSC_USE_COMPLEX) 4686 for (k=0; k<a->nz; k++) { 4687 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4688 *flg = PETSC_FALSE; 4689 PetscFunctionReturn(0); 4690 } 4691 } 4692 #else 4693 ierr = PetscArraycmp(a->a,b->a,a->nz,flg);CHKERRQ(ierr); 4694 #endif 4695 PetscFunctionReturn(0); 4696 } 4697 4698 /*@ 4699 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4700 provided by the user. 4701 4702 Collective 4703 4704 Input Parameters: 4705 + comm - must be an MPI communicator of size 1 4706 . m - number of rows 4707 . n - number of columns 4708 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 4709 . j - column indices 4710 - a - matrix values 4711 4712 Output Parameter: 4713 . mat - the matrix 4714 4715 Level: intermediate 4716 4717 Notes: 4718 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4719 once the matrix is destroyed and not before 4720 4721 You cannot set new nonzero locations into this matrix, that will generate an error. 4722 4723 The i and j indices are 0 based 4724 4725 The format which is used for the sparse matrix input, is equivalent to a 4726 row-major ordering.. i.e for the following matrix, the input data expected is 4727 as shown 4728 4729 $ 1 0 0 4730 $ 2 0 3 4731 $ 4 5 6 4732 $ 4733 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4734 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4735 $ v = {1,2,3,4,5,6} [size = 6] 4736 4737 4738 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4739 4740 @*/ 4741 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4742 { 4743 PetscErrorCode ierr; 4744 PetscInt ii; 4745 Mat_SeqAIJ *aij; 4746 PetscInt jj; 4747 4748 PetscFunctionBegin; 4749 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4750 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4751 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4752 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4753 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4754 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 4755 aij = (Mat_SeqAIJ*)(*mat)->data; 4756 ierr = PetscMalloc1(m,&aij->imax);CHKERRQ(ierr); 4757 ierr = PetscMalloc1(m,&aij->ilen);CHKERRQ(ierr); 4758 4759 aij->i = i; 4760 aij->j = j; 4761 aij->a = a; 4762 aij->singlemalloc = PETSC_FALSE; 4763 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4764 aij->free_a = PETSC_FALSE; 4765 aij->free_ij = PETSC_FALSE; 4766 4767 for (ii=0; ii<m; ii++) { 4768 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4769 if (PetscDefined(USE_DEBUG)) { 4770 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]); 4771 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4772 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); 4773 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); 4774 } 4775 } 4776 } 4777 if (PetscDefined(USE_DEBUG)) { 4778 for (ii=0; ii<aij->i[m]; ii++) { 4779 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4780 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]); 4781 } 4782 } 4783 4784 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4785 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4786 PetscFunctionReturn(0); 4787 } 4788 /*@C 4789 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4790 provided by the user. 4791 4792 Collective 4793 4794 Input Parameters: 4795 + comm - must be an MPI communicator of size 1 4796 . m - number of rows 4797 . n - number of columns 4798 . i - row indices 4799 . j - column indices 4800 . a - matrix values 4801 . nz - number of nonzeros 4802 - idx - 0 or 1 based 4803 4804 Output Parameter: 4805 . mat - the matrix 4806 4807 Level: intermediate 4808 4809 Notes: 4810 The i and j indices are 0 based 4811 4812 The format which is used for the sparse matrix input, is equivalent to a 4813 row-major ordering.. i.e for the following matrix, the input data expected is 4814 as shown: 4815 4816 1 0 0 4817 2 0 3 4818 4 5 6 4819 4820 i = {0,1,1,2,2,2} 4821 j = {0,0,2,0,1,2} 4822 v = {1,2,3,4,5,6} 4823 4824 4825 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4826 4827 @*/ 4828 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4829 { 4830 PetscErrorCode ierr; 4831 PetscInt ii, *nnz, one = 1,row,col; 4832 4833 4834 PetscFunctionBegin; 4835 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4836 for (ii = 0; ii < nz; ii++) { 4837 nnz[i[ii] - !!idx] += 1; 4838 } 4839 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4840 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4841 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4842 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4843 for (ii = 0; ii < nz; ii++) { 4844 if (idx) { 4845 row = i[ii] - 1; 4846 col = j[ii] - 1; 4847 } else { 4848 row = i[ii]; 4849 col = j[ii]; 4850 } 4851 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4852 } 4853 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4854 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4855 ierr = PetscFree(nnz);CHKERRQ(ierr); 4856 PetscFunctionReturn(0); 4857 } 4858 4859 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4860 { 4861 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4862 PetscErrorCode ierr; 4863 4864 PetscFunctionBegin; 4865 a->idiagvalid = PETSC_FALSE; 4866 a->ibdiagvalid = PETSC_FALSE; 4867 4868 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4869 PetscFunctionReturn(0); 4870 } 4871 4872 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4873 { 4874 PetscErrorCode ierr; 4875 PetscMPIInt size; 4876 4877 PetscFunctionBegin; 4878 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4879 if (size == 1) { 4880 if (scall == MAT_INITIAL_MATRIX) { 4881 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 4882 } else { 4883 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4884 } 4885 } else { 4886 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4887 } 4888 PetscFunctionReturn(0); 4889 } 4890 4891 /* 4892 Permute A into C's *local* index space using rowemb,colemb. 4893 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4894 of [0,m), colemb is in [0,n). 4895 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4896 */ 4897 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4898 { 4899 /* If making this function public, change the error returned in this function away from _PLIB. */ 4900 PetscErrorCode ierr; 4901 Mat_SeqAIJ *Baij; 4902 PetscBool seqaij; 4903 PetscInt m,n,*nz,i,j,count; 4904 PetscScalar v; 4905 const PetscInt *rowindices,*colindices; 4906 4907 PetscFunctionBegin; 4908 if (!B) PetscFunctionReturn(0); 4909 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4910 ierr = PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4911 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4912 if (rowemb) { 4913 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4914 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); 4915 } else { 4916 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4917 } 4918 if (colemb) { 4919 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4920 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); 4921 } else { 4922 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4923 } 4924 4925 Baij = (Mat_SeqAIJ*)(B->data); 4926 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4927 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4928 for (i=0; i<B->rmap->n; i++) { 4929 nz[i] = Baij->i[i+1] - Baij->i[i]; 4930 } 4931 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4932 ierr = PetscFree(nz);CHKERRQ(ierr); 4933 } 4934 if (pattern == SUBSET_NONZERO_PATTERN) { 4935 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4936 } 4937 count = 0; 4938 rowindices = NULL; 4939 colindices = NULL; 4940 if (rowemb) { 4941 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4942 } 4943 if (colemb) { 4944 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4945 } 4946 for (i=0; i<B->rmap->n; i++) { 4947 PetscInt row; 4948 row = i; 4949 if (rowindices) row = rowindices[i]; 4950 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4951 PetscInt col; 4952 col = Baij->j[count]; 4953 if (colindices) col = colindices[col]; 4954 v = Baij->a[count]; 4955 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4956 ++count; 4957 } 4958 } 4959 /* FIXME: set C's nonzerostate correctly. */ 4960 /* Assembly for C is necessary. */ 4961 C->preallocated = PETSC_TRUE; 4962 C->assembled = PETSC_TRUE; 4963 C->was_assembled = PETSC_FALSE; 4964 PetscFunctionReturn(0); 4965 } 4966 4967 PetscFunctionList MatSeqAIJList = NULL; 4968 4969 /*@C 4970 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 4971 4972 Collective on Mat 4973 4974 Input Parameters: 4975 + mat - the matrix object 4976 - matype - matrix type 4977 4978 Options Database Key: 4979 . -mat_seqai_type <method> - for example seqaijcrl 4980 4981 4982 Level: intermediate 4983 4984 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 4985 @*/ 4986 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 4987 { 4988 PetscErrorCode ierr,(*r)(Mat,MatType,MatReuse,Mat*); 4989 PetscBool sametype; 4990 4991 PetscFunctionBegin; 4992 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 4993 ierr = PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);CHKERRQ(ierr); 4994 if (sametype) PetscFunctionReturn(0); 4995 4996 ierr = PetscFunctionListFind(MatSeqAIJList,matype,&r);CHKERRQ(ierr); 4997 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 4998 ierr = (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);CHKERRQ(ierr); 4999 PetscFunctionReturn(0); 5000 } 5001 5002 5003 /*@C 5004 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 5005 5006 Not Collective 5007 5008 Input Parameters: 5009 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 5010 - function - routine to convert to subtype 5011 5012 Notes: 5013 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 5014 5015 5016 Then, your matrix can be chosen with the procedural interface at runtime via the option 5017 $ -mat_seqaij_type my_mat 5018 5019 Level: advanced 5020 5021 .seealso: MatSeqAIJRegisterAll() 5022 5023 5024 Level: advanced 5025 @*/ 5026 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 5027 { 5028 PetscErrorCode ierr; 5029 5030 PetscFunctionBegin; 5031 ierr = MatInitializePackage();CHKERRQ(ierr); 5032 ierr = PetscFunctionListAdd(&MatSeqAIJList,sname,function);CHKERRQ(ierr); 5033 PetscFunctionReturn(0); 5034 } 5035 5036 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5037 5038 /*@C 5039 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 5040 5041 Not Collective 5042 5043 Level: advanced 5044 5045 Developers Note: CUSP and CUSPARSE do not yet support the MatConvert_SeqAIJ..() paradigm and thus cannot be registered here 5046 5047 .seealso: MatRegisterAll(), MatSeqAIJRegister() 5048 @*/ 5049 PetscErrorCode MatSeqAIJRegisterAll(void) 5050 { 5051 PetscErrorCode ierr; 5052 5053 PetscFunctionBegin; 5054 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5055 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5056 5057 ierr = MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 5058 ierr = MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 5059 ierr = MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL);CHKERRQ(ierr); 5060 #if defined(PETSC_HAVE_MKL_SPARSE) 5061 ierr = MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 5062 #endif 5063 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5064 ierr = MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr); 5065 #endif 5066 PetscFunctionReturn(0); 5067 } 5068 5069 /* 5070 Special version for direct calls from Fortran 5071 */ 5072 #include <petsc/private/fortranimpl.h> 5073 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5074 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5075 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5076 #define matsetvaluesseqaij_ matsetvaluesseqaij 5077 #endif 5078 5079 /* Change these macros so can be used in void function */ 5080 #undef CHKERRQ 5081 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 5082 #undef SETERRQ2 5083 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 5084 #undef SETERRQ3 5085 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 5086 5087 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 5088 { 5089 Mat A = *AA; 5090 PetscInt m = *mm, n = *nn; 5091 InsertMode is = *isis; 5092 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 5093 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 5094 PetscInt *imax,*ai,*ailen; 5095 PetscErrorCode ierr; 5096 PetscInt *aj,nonew = a->nonew,lastcol = -1; 5097 MatScalar *ap,value,*aa; 5098 PetscBool ignorezeroentries = a->ignorezeroentries; 5099 PetscBool roworiented = a->roworiented; 5100 5101 PetscFunctionBegin; 5102 MatCheckPreallocated(A,1); 5103 imax = a->imax; 5104 ai = a->i; 5105 ailen = a->ilen; 5106 aj = a->j; 5107 aa = a->a; 5108 5109 for (k=0; k<m; k++) { /* loop over added rows */ 5110 row = im[k]; 5111 if (row < 0) continue; 5112 if (PetscUnlikelyDebug(row >= A->rmap->n)) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 5113 rp = aj + ai[row]; ap = aa + ai[row]; 5114 rmax = imax[row]; nrow = ailen[row]; 5115 low = 0; 5116 high = nrow; 5117 for (l=0; l<n; l++) { /* loop over added columns */ 5118 if (in[l] < 0) continue; 5119 if (PetscUnlikelyDebug(in[l] >= A->cmap->n)) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 5120 col = in[l]; 5121 if (roworiented) value = v[l + k*n]; 5122 else value = v[k + l*m]; 5123 5124 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5125 5126 if (col <= lastcol) low = 0; 5127 else high = nrow; 5128 lastcol = col; 5129 while (high-low > 5) { 5130 t = (low+high)/2; 5131 if (rp[t] > col) high = t; 5132 else low = t; 5133 } 5134 for (i=low; i<high; i++) { 5135 if (rp[i] > col) break; 5136 if (rp[i] == col) { 5137 if (is == ADD_VALUES) ap[i] += value; 5138 else ap[i] = value; 5139 goto noinsert; 5140 } 5141 } 5142 if (value == 0.0 && ignorezeroentries) goto noinsert; 5143 if (nonew == 1) goto noinsert; 5144 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 5145 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 5146 N = nrow++ - 1; a->nz++; high++; 5147 /* shift up all the later entries in this row */ 5148 for (ii=N; ii>=i; ii--) { 5149 rp[ii+1] = rp[ii]; 5150 ap[ii+1] = ap[ii]; 5151 } 5152 rp[i] = col; 5153 ap[i] = value; 5154 A->nonzerostate++; 5155 noinsert:; 5156 low = i + 1; 5157 } 5158 ailen[row] = nrow; 5159 } 5160 PetscFunctionReturnVoid(); 5161 } 5162