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