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