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