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