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