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+1,&collengths);CHKERRQ(ierr); 261 ierr = PetscMalloc1(n+1,&cia);CHKERRQ(ierr); 262 ierr = PetscMalloc1(nz+1,&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+1,&collengths);CHKERRQ(ierr); 318 ierr = PetscMalloc1(n+1,&cia);CHKERRQ(ierr); 319 ierr = PetscMalloc1(nz+1,&cja);CHKERRQ(ierr); 320 ierr = PetscMalloc1(nz+1,&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 2207 b = (Mat_SeqAIJ*)((*B)->data); 2208 b->free_a = PETSC_FALSE; 2209 b->free_ij = PETSC_TRUE; 2210 b->nonew = 0; 2211 PetscFunctionReturn(0); 2212 } 2213 2214 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2215 { 2216 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2217 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2218 MatScalar *va,*vb; 2219 PetscErrorCode ierr; 2220 PetscInt ma,na,mb,nb, i; 2221 2222 PetscFunctionBegin; 2223 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2224 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2225 if (ma!=nb || na!=mb) { 2226 *f = PETSC_FALSE; 2227 PetscFunctionReturn(0); 2228 } 2229 aii = aij->i; bii = bij->i; 2230 adx = aij->j; bdx = bij->j; 2231 va = aij->a; vb = bij->a; 2232 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2233 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2234 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2235 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2236 2237 *f = PETSC_TRUE; 2238 for (i=0; i<ma; i++) { 2239 while (aptr[i]<aii[i+1]) { 2240 PetscInt idc,idr; 2241 PetscScalar vc,vr; 2242 /* column/row index/value */ 2243 idc = adx[aptr[i]]; 2244 idr = bdx[bptr[idc]]; 2245 vc = va[aptr[i]]; 2246 vr = vb[bptr[idc]]; 2247 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2248 *f = PETSC_FALSE; 2249 goto done; 2250 } else { 2251 aptr[i]++; 2252 if (B || i!=idc) bptr[idc]++; 2253 } 2254 } 2255 } 2256 done: 2257 ierr = PetscFree(aptr);CHKERRQ(ierr); 2258 ierr = PetscFree(bptr);CHKERRQ(ierr); 2259 PetscFunctionReturn(0); 2260 } 2261 2262 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2263 { 2264 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2265 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2266 MatScalar *va,*vb; 2267 PetscErrorCode ierr; 2268 PetscInt ma,na,mb,nb, i; 2269 2270 PetscFunctionBegin; 2271 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2272 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2273 if (ma!=nb || na!=mb) { 2274 *f = PETSC_FALSE; 2275 PetscFunctionReturn(0); 2276 } 2277 aii = aij->i; bii = bij->i; 2278 adx = aij->j; bdx = bij->j; 2279 va = aij->a; vb = bij->a; 2280 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2281 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2282 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2283 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2284 2285 *f = PETSC_TRUE; 2286 for (i=0; i<ma; i++) { 2287 while (aptr[i]<aii[i+1]) { 2288 PetscInt idc,idr; 2289 PetscScalar vc,vr; 2290 /* column/row index/value */ 2291 idc = adx[aptr[i]]; 2292 idr = bdx[bptr[idc]]; 2293 vc = va[aptr[i]]; 2294 vr = vb[bptr[idc]]; 2295 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2296 *f = PETSC_FALSE; 2297 goto done; 2298 } else { 2299 aptr[i]++; 2300 if (B || i!=idc) bptr[idc]++; 2301 } 2302 } 2303 } 2304 done: 2305 ierr = PetscFree(aptr);CHKERRQ(ierr); 2306 ierr = PetscFree(bptr);CHKERRQ(ierr); 2307 PetscFunctionReturn(0); 2308 } 2309 2310 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2311 { 2312 PetscErrorCode ierr; 2313 2314 PetscFunctionBegin; 2315 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2316 PetscFunctionReturn(0); 2317 } 2318 2319 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2320 { 2321 PetscErrorCode ierr; 2322 2323 PetscFunctionBegin; 2324 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2325 PetscFunctionReturn(0); 2326 } 2327 2328 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2329 { 2330 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2331 const PetscScalar *l,*r; 2332 PetscScalar x; 2333 MatScalar *v; 2334 PetscErrorCode ierr; 2335 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2336 const PetscInt *jj; 2337 2338 PetscFunctionBegin; 2339 if (ll) { 2340 /* The local size is used so that VecMPI can be passed to this routine 2341 by MatDiagonalScale_MPIAIJ */ 2342 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2343 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2344 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 2345 v = a->a; 2346 for (i=0; i<m; i++) { 2347 x = l[i]; 2348 M = a->i[i+1] - a->i[i]; 2349 for (j=0; j<M; j++) (*v++) *= x; 2350 } 2351 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 2352 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2353 } 2354 if (rr) { 2355 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2356 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2357 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 2358 v = a->a; jj = a->j; 2359 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2360 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 2361 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2362 } 2363 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2364 PetscFunctionReturn(0); 2365 } 2366 2367 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2368 { 2369 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2370 PetscErrorCode ierr; 2371 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2372 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2373 const PetscInt *irow,*icol; 2374 PetscInt nrows,ncols; 2375 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2376 MatScalar *a_new,*mat_a; 2377 Mat C; 2378 PetscBool stride; 2379 2380 PetscFunctionBegin; 2381 2382 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2383 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2384 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2385 2386 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2387 if (stride) { 2388 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2389 } else { 2390 first = 0; 2391 step = 0; 2392 } 2393 if (stride && step == 1) { 2394 /* special case of contiguous rows */ 2395 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2396 /* loop over new rows determining lens and starting points */ 2397 for (i=0; i<nrows; i++) { 2398 kstart = ai[irow[i]]; 2399 kend = kstart + ailen[irow[i]]; 2400 starts[i] = kstart; 2401 for (k=kstart; k<kend; k++) { 2402 if (aj[k] >= first) { 2403 starts[i] = k; 2404 break; 2405 } 2406 } 2407 sum = 0; 2408 while (k < kend) { 2409 if (aj[k++] >= first+ncols) break; 2410 sum++; 2411 } 2412 lens[i] = sum; 2413 } 2414 /* create submatrix */ 2415 if (scall == MAT_REUSE_MATRIX) { 2416 PetscInt n_cols,n_rows; 2417 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2418 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2419 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2420 C = *B; 2421 } else { 2422 PetscInt rbs,cbs; 2423 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2424 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2425 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2426 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2427 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2428 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2429 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2430 } 2431 c = (Mat_SeqAIJ*)C->data; 2432 2433 /* loop over rows inserting into submatrix */ 2434 a_new = c->a; 2435 j_new = c->j; 2436 i_new = c->i; 2437 2438 for (i=0; i<nrows; i++) { 2439 ii = starts[i]; 2440 lensi = lens[i]; 2441 for (k=0; k<lensi; k++) { 2442 *j_new++ = aj[ii+k] - first; 2443 } 2444 ierr = PetscArraycpy(a_new,a->a + starts[i],lensi);CHKERRQ(ierr); 2445 a_new += lensi; 2446 i_new[i+1] = i_new[i] + lensi; 2447 c->ilen[i] = lensi; 2448 } 2449 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2450 } else { 2451 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2452 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2453 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2454 for (i=0; i<ncols; i++) { 2455 #if defined(PETSC_USE_DEBUG) 2456 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); 2457 #endif 2458 smap[icol[i]] = i+1; 2459 } 2460 2461 /* determine lens of each row */ 2462 for (i=0; i<nrows; i++) { 2463 kstart = ai[irow[i]]; 2464 kend = kstart + a->ilen[irow[i]]; 2465 lens[i] = 0; 2466 for (k=kstart; k<kend; k++) { 2467 if (smap[aj[k]]) { 2468 lens[i]++; 2469 } 2470 } 2471 } 2472 /* Create and fill new matrix */ 2473 if (scall == MAT_REUSE_MATRIX) { 2474 PetscBool equal; 2475 2476 c = (Mat_SeqAIJ*)((*B)->data); 2477 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2478 ierr = PetscArraycmp(c->ilen,lens,(*B)->rmap->n,&equal);CHKERRQ(ierr); 2479 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2480 ierr = PetscArrayzero(c->ilen,(*B)->rmap->n);CHKERRQ(ierr); 2481 C = *B; 2482 } else { 2483 PetscInt rbs,cbs; 2484 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2485 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2486 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2487 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2488 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2489 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2490 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2491 } 2492 c = (Mat_SeqAIJ*)(C->data); 2493 for (i=0; i<nrows; i++) { 2494 row = irow[i]; 2495 kstart = ai[row]; 2496 kend = kstart + a->ilen[row]; 2497 mat_i = c->i[i]; 2498 mat_j = c->j + mat_i; 2499 mat_a = c->a + mat_i; 2500 mat_ilen = c->ilen + i; 2501 for (k=kstart; k<kend; k++) { 2502 if ((tcol=smap[a->j[k]])) { 2503 *mat_j++ = tcol - 1; 2504 *mat_a++ = a->a[k]; 2505 (*mat_ilen)++; 2506 2507 } 2508 } 2509 } 2510 /* Free work space */ 2511 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2512 ierr = PetscFree(smap);CHKERRQ(ierr); 2513 ierr = PetscFree(lens);CHKERRQ(ierr); 2514 /* sort */ 2515 for (i = 0; i < nrows; i++) { 2516 PetscInt ilen; 2517 2518 mat_i = c->i[i]; 2519 mat_j = c->j + mat_i; 2520 mat_a = c->a + mat_i; 2521 ilen = c->ilen[i]; 2522 ierr = PetscSortIntWithScalarArray(ilen,mat_j,mat_a);CHKERRQ(ierr); 2523 } 2524 } 2525 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2526 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2527 2528 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2529 *B = C; 2530 PetscFunctionReturn(0); 2531 } 2532 2533 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2534 { 2535 PetscErrorCode ierr; 2536 Mat B; 2537 2538 PetscFunctionBegin; 2539 if (scall == MAT_INITIAL_MATRIX) { 2540 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2541 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2542 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2543 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2544 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2545 *subMat = B; 2546 } else { 2547 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2548 } 2549 PetscFunctionReturn(0); 2550 } 2551 2552 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2553 { 2554 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2555 PetscErrorCode ierr; 2556 Mat outA; 2557 PetscBool row_identity,col_identity; 2558 2559 PetscFunctionBegin; 2560 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2561 2562 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2563 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2564 2565 outA = inA; 2566 outA->factortype = MAT_FACTOR_LU; 2567 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 2568 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 2569 2570 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2571 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2572 2573 a->row = row; 2574 2575 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2576 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2577 2578 a->col = col; 2579 2580 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2581 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2582 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2583 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2584 2585 if (!a->solve_work) { /* this matrix may have been factored before */ 2586 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2587 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2588 } 2589 2590 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2591 if (row_identity && col_identity) { 2592 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2593 } else { 2594 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2595 } 2596 PetscFunctionReturn(0); 2597 } 2598 2599 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2600 { 2601 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2602 PetscScalar oalpha = alpha; 2603 PetscErrorCode ierr; 2604 PetscBLASInt one = 1,bnz; 2605 2606 PetscFunctionBegin; 2607 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2608 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2609 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2610 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2611 PetscFunctionReturn(0); 2612 } 2613 2614 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2615 { 2616 PetscErrorCode ierr; 2617 PetscInt i; 2618 2619 PetscFunctionBegin; 2620 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2621 ierr = PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr);CHKERRQ(ierr); 2622 2623 for (i=0; i<submatj->nrqr; ++i) { 2624 ierr = PetscFree(submatj->sbuf2[i]);CHKERRQ(ierr); 2625 } 2626 ierr = PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1);CHKERRQ(ierr); 2627 2628 if (submatj->rbuf1) { 2629 ierr = PetscFree(submatj->rbuf1[0]);CHKERRQ(ierr); 2630 ierr = PetscFree(submatj->rbuf1);CHKERRQ(ierr); 2631 } 2632 2633 for (i=0; i<submatj->nrqs; ++i) { 2634 ierr = PetscFree(submatj->rbuf3[i]);CHKERRQ(ierr); 2635 } 2636 ierr = PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3);CHKERRQ(ierr); 2637 ierr = PetscFree(submatj->pa);CHKERRQ(ierr); 2638 } 2639 2640 #if defined(PETSC_USE_CTABLE) 2641 ierr = PetscTableDestroy((PetscTable*)&submatj->rmap);CHKERRQ(ierr); 2642 if (submatj->cmap_loc) {ierr = PetscFree(submatj->cmap_loc);CHKERRQ(ierr);} 2643 ierr = PetscFree(submatj->rmap_loc);CHKERRQ(ierr); 2644 #else 2645 ierr = PetscFree(submatj->rmap);CHKERRQ(ierr); 2646 #endif 2647 2648 if (!submatj->allcolumns) { 2649 #if defined(PETSC_USE_CTABLE) 2650 ierr = PetscTableDestroy((PetscTable*)&submatj->cmap);CHKERRQ(ierr); 2651 #else 2652 ierr = PetscFree(submatj->cmap);CHKERRQ(ierr); 2653 #endif 2654 } 2655 ierr = PetscFree(submatj->row2proc);CHKERRQ(ierr); 2656 2657 ierr = PetscFree(submatj);CHKERRQ(ierr); 2658 PetscFunctionReturn(0); 2659 } 2660 2661 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2662 { 2663 PetscErrorCode ierr; 2664 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2665 Mat_SubSppt *submatj = c->submatis1; 2666 2667 PetscFunctionBegin; 2668 ierr = (*submatj->destroy)(C);CHKERRQ(ierr); 2669 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2670 PetscFunctionReturn(0); 2671 } 2672 2673 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2674 { 2675 PetscErrorCode ierr; 2676 PetscInt i; 2677 Mat C; 2678 Mat_SeqAIJ *c; 2679 Mat_SubSppt *submatj; 2680 2681 PetscFunctionBegin; 2682 for (i=0; i<n; i++) { 2683 C = (*mat)[i]; 2684 c = (Mat_SeqAIJ*)C->data; 2685 submatj = c->submatis1; 2686 if (submatj) { 2687 if (--((PetscObject)C)->refct <= 0) { 2688 ierr = (*submatj->destroy)(C);CHKERRQ(ierr); 2689 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2690 ierr = PetscFree(C->defaultvectype);CHKERRQ(ierr); 2691 ierr = PetscLayoutDestroy(&C->rmap);CHKERRQ(ierr); 2692 ierr = PetscLayoutDestroy(&C->cmap);CHKERRQ(ierr); 2693 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 2694 } 2695 } else { 2696 ierr = MatDestroy(&C);CHKERRQ(ierr); 2697 } 2698 } 2699 2700 /* Destroy Dummy submatrices created for reuse */ 2701 ierr = MatDestroySubMatrices_Dummy(n,mat);CHKERRQ(ierr); 2702 2703 ierr = PetscFree(*mat);CHKERRQ(ierr); 2704 PetscFunctionReturn(0); 2705 } 2706 2707 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2708 { 2709 PetscErrorCode ierr; 2710 PetscInt i; 2711 2712 PetscFunctionBegin; 2713 if (scall == MAT_INITIAL_MATRIX) { 2714 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2715 } 2716 2717 for (i=0; i<n; i++) { 2718 ierr = MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2719 } 2720 PetscFunctionReturn(0); 2721 } 2722 2723 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2724 { 2725 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2726 PetscErrorCode ierr; 2727 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2728 const PetscInt *idx; 2729 PetscInt start,end,*ai,*aj; 2730 PetscBT table; 2731 2732 PetscFunctionBegin; 2733 m = A->rmap->n; 2734 ai = a->i; 2735 aj = a->j; 2736 2737 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2738 2739 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2740 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2741 2742 for (i=0; i<is_max; i++) { 2743 /* Initialize the two local arrays */ 2744 isz = 0; 2745 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2746 2747 /* Extract the indices, assume there can be duplicate entries */ 2748 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2749 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2750 2751 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2752 for (j=0; j<n; ++j) { 2753 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2754 } 2755 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2756 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2757 2758 k = 0; 2759 for (j=0; j<ov; j++) { /* for each overlap */ 2760 n = isz; 2761 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2762 row = nidx[k]; 2763 start = ai[row]; 2764 end = ai[row+1]; 2765 for (l = start; l<end; l++) { 2766 val = aj[l]; 2767 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2768 } 2769 } 2770 } 2771 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2772 } 2773 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2774 ierr = PetscFree(nidx);CHKERRQ(ierr); 2775 PetscFunctionReturn(0); 2776 } 2777 2778 /* -------------------------------------------------------------- */ 2779 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2780 { 2781 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2782 PetscErrorCode ierr; 2783 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2784 const PetscInt *row,*col; 2785 PetscInt *cnew,j,*lens; 2786 IS icolp,irowp; 2787 PetscInt *cwork = NULL; 2788 PetscScalar *vwork = NULL; 2789 2790 PetscFunctionBegin; 2791 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2792 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2793 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2794 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2795 2796 /* determine lengths of permuted rows */ 2797 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2798 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2799 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2800 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2801 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2802 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2803 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2804 ierr = PetscFree(lens);CHKERRQ(ierr); 2805 2806 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2807 for (i=0; i<m; i++) { 2808 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2809 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2810 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2811 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2812 } 2813 ierr = PetscFree(cnew);CHKERRQ(ierr); 2814 2815 (*B)->assembled = PETSC_FALSE; 2816 2817 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2818 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2819 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2820 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2821 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2822 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2823 PetscFunctionReturn(0); 2824 } 2825 2826 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2827 { 2828 PetscErrorCode ierr; 2829 2830 PetscFunctionBegin; 2831 /* If the two matrices have the same copy implementation, use fast copy. */ 2832 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2833 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2834 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2835 2836 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"); 2837 ierr = PetscArraycpy(b->a,a->a,a->i[A->rmap->n]);CHKERRQ(ierr); 2838 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2839 } else { 2840 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2841 } 2842 PetscFunctionReturn(0); 2843 } 2844 2845 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2846 { 2847 PetscErrorCode ierr; 2848 2849 PetscFunctionBegin; 2850 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2851 PetscFunctionReturn(0); 2852 } 2853 2854 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2855 { 2856 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2857 2858 PetscFunctionBegin; 2859 *array = a->a; 2860 PetscFunctionReturn(0); 2861 } 2862 2863 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2864 { 2865 PetscFunctionBegin; 2866 PetscFunctionReturn(0); 2867 } 2868 2869 /* 2870 Computes the number of nonzeros per row needed for preallocation when X and Y 2871 have different nonzero structure. 2872 */ 2873 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2874 { 2875 PetscInt i,j,k,nzx,nzy; 2876 2877 PetscFunctionBegin; 2878 /* Set the number of nonzeros in the new matrix */ 2879 for (i=0; i<m; i++) { 2880 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2881 nzx = xi[i+1] - xi[i]; 2882 nzy = yi[i+1] - yi[i]; 2883 nnz[i] = 0; 2884 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2885 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2886 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2887 nnz[i]++; 2888 } 2889 for (; k<nzy; k++) nnz[i]++; 2890 } 2891 PetscFunctionReturn(0); 2892 } 2893 2894 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2895 { 2896 PetscInt m = Y->rmap->N; 2897 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2898 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2899 PetscErrorCode ierr; 2900 2901 PetscFunctionBegin; 2902 /* Set the number of nonzeros in the new matrix */ 2903 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2904 PetscFunctionReturn(0); 2905 } 2906 2907 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2908 { 2909 PetscErrorCode ierr; 2910 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2911 PetscBLASInt one=1,bnz; 2912 2913 PetscFunctionBegin; 2914 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2915 if (str == SAME_NONZERO_PATTERN) { 2916 PetscScalar alpha = a; 2917 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2918 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2919 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2920 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2921 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2922 } else { 2923 Mat B; 2924 PetscInt *nnz; 2925 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2926 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2927 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2928 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2929 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2930 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2931 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2932 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2933 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2934 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2935 ierr = PetscFree(nnz);CHKERRQ(ierr); 2936 } 2937 PetscFunctionReturn(0); 2938 } 2939 2940 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2941 { 2942 #if defined(PETSC_USE_COMPLEX) 2943 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2944 PetscInt i,nz; 2945 PetscScalar *a; 2946 2947 PetscFunctionBegin; 2948 nz = aij->nz; 2949 a = aij->a; 2950 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2951 #else 2952 PetscFunctionBegin; 2953 #endif 2954 PetscFunctionReturn(0); 2955 } 2956 2957 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2958 { 2959 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2960 PetscErrorCode ierr; 2961 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2962 PetscReal atmp; 2963 PetscScalar *x; 2964 MatScalar *aa; 2965 2966 PetscFunctionBegin; 2967 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2968 aa = a->a; 2969 ai = a->i; 2970 aj = a->j; 2971 2972 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2973 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2974 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2975 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2976 for (i=0; i<m; i++) { 2977 ncols = ai[1] - ai[0]; ai++; 2978 x[i] = 0.0; 2979 for (j=0; j<ncols; j++) { 2980 atmp = PetscAbsScalar(*aa); 2981 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2982 aa++; aj++; 2983 } 2984 } 2985 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2986 PetscFunctionReturn(0); 2987 } 2988 2989 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2990 { 2991 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2992 PetscErrorCode ierr; 2993 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2994 PetscScalar *x; 2995 MatScalar *aa; 2996 2997 PetscFunctionBegin; 2998 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2999 aa = a->a; 3000 ai = a->i; 3001 aj = a->j; 3002 3003 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3004 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3005 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3006 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3007 for (i=0; i<m; i++) { 3008 ncols = ai[1] - ai[0]; ai++; 3009 if (ncols == A->cmap->n) { /* row is dense */ 3010 x[i] = *aa; if (idx) idx[i] = 0; 3011 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3012 x[i] = 0.0; 3013 if (idx) { 3014 idx[i] = 0; /* in case ncols is zero */ 3015 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 3016 if (aj[j] > j) { 3017 idx[i] = j; 3018 break; 3019 } 3020 } 3021 } 3022 } 3023 for (j=0; j<ncols; j++) { 3024 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3025 aa++; aj++; 3026 } 3027 } 3028 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3029 PetscFunctionReturn(0); 3030 } 3031 3032 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3033 { 3034 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3035 PetscErrorCode ierr; 3036 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3037 PetscReal atmp; 3038 PetscScalar *x; 3039 MatScalar *aa; 3040 3041 PetscFunctionBegin; 3042 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3043 aa = a->a; 3044 ai = a->i; 3045 aj = a->j; 3046 3047 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3048 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3049 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3050 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); 3051 for (i=0; i<m; i++) { 3052 ncols = ai[1] - ai[0]; ai++; 3053 if (ncols) { 3054 /* Get first nonzero */ 3055 for (j = 0; j < ncols; j++) { 3056 atmp = PetscAbsScalar(aa[j]); 3057 if (atmp > 1.0e-12) { 3058 x[i] = atmp; 3059 if (idx) idx[i] = aj[j]; 3060 break; 3061 } 3062 } 3063 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 3064 } else { 3065 x[i] = 0.0; if (idx) idx[i] = 0; 3066 } 3067 for (j = 0; j < ncols; j++) { 3068 atmp = PetscAbsScalar(*aa); 3069 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3070 aa++; aj++; 3071 } 3072 } 3073 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3074 PetscFunctionReturn(0); 3075 } 3076 3077 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3078 { 3079 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3080 PetscErrorCode ierr; 3081 PetscInt i,j,m = A->rmap->n,ncols,n; 3082 const PetscInt *ai,*aj; 3083 PetscScalar *x; 3084 const MatScalar *aa; 3085 3086 PetscFunctionBegin; 3087 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3088 aa = a->a; 3089 ai = a->i; 3090 aj = a->j; 3091 3092 ierr = VecSet(v,0.0);CHKERRQ(ierr); 3093 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 3094 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 3095 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3096 for (i=0; i<m; i++) { 3097 ncols = ai[1] - ai[0]; ai++; 3098 if (ncols == A->cmap->n) { /* row is dense */ 3099 x[i] = *aa; if (idx) idx[i] = 0; 3100 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3101 x[i] = 0.0; 3102 if (idx) { /* find first implicit 0.0 in the row */ 3103 idx[i] = 0; /* in case ncols is zero */ 3104 for (j=0; j<ncols; j++) { 3105 if (aj[j] > j) { 3106 idx[i] = j; 3107 break; 3108 } 3109 } 3110 } 3111 } 3112 for (j=0; j<ncols; j++) { 3113 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3114 aa++; aj++; 3115 } 3116 } 3117 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 3118 PetscFunctionReturn(0); 3119 } 3120 3121 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3122 { 3123 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3124 PetscErrorCode ierr; 3125 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3126 MatScalar *diag,work[25],*v_work; 3127 const PetscReal shift = 0.0; 3128 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3129 3130 PetscFunctionBegin; 3131 allowzeropivot = PetscNot(A->erroriffailure); 3132 if (a->ibdiagvalid) { 3133 if (values) *values = a->ibdiag; 3134 PetscFunctionReturn(0); 3135 } 3136 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 3137 if (!a->ibdiag) { 3138 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 3139 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 3140 } 3141 diag = a->ibdiag; 3142 if (values) *values = a->ibdiag; 3143 /* factor and invert each block */ 3144 switch (bs) { 3145 case 1: 3146 for (i=0; i<mbs; i++) { 3147 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3148 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3149 if (allowzeropivot) { 3150 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3151 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3152 A->factorerror_zeropivot_row = i; 3153 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 3154 } 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); 3155 } 3156 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3157 } 3158 break; 3159 case 2: 3160 for (i=0; i<mbs; i++) { 3161 ij[0] = 2*i; ij[1] = 2*i + 1; 3162 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3163 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3164 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3165 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3166 diag += 4; 3167 } 3168 break; 3169 case 3: 3170 for (i=0; i<mbs; i++) { 3171 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3172 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3173 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3174 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3175 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3176 diag += 9; 3177 } 3178 break; 3179 case 4: 3180 for (i=0; i<mbs; i++) { 3181 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3182 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3183 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3184 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3185 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3186 diag += 16; 3187 } 3188 break; 3189 case 5: 3190 for (i=0; i<mbs; i++) { 3191 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3192 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3193 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3194 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3195 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3196 diag += 25; 3197 } 3198 break; 3199 case 6: 3200 for (i=0; i<mbs; i++) { 3201 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; 3202 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3203 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3204 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3205 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3206 diag += 36; 3207 } 3208 break; 3209 case 7: 3210 for (i=0; i<mbs; i++) { 3211 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; 3212 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3213 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3214 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3215 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3216 diag += 49; 3217 } 3218 break; 3219 default: 3220 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3221 for (i=0; i<mbs; i++) { 3222 for (j=0; j<bs; j++) { 3223 IJ[j] = bs*i + j; 3224 } 3225 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3226 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3227 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3228 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3229 diag += bs2; 3230 } 3231 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3232 } 3233 a->ibdiagvalid = PETSC_TRUE; 3234 PetscFunctionReturn(0); 3235 } 3236 3237 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3238 { 3239 PetscErrorCode ierr; 3240 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3241 PetscScalar a; 3242 PetscInt m,n,i,j,col; 3243 3244 PetscFunctionBegin; 3245 if (!x->assembled) { 3246 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3247 for (i=0; i<m; i++) { 3248 for (j=0; j<aij->imax[i]; j++) { 3249 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3250 col = (PetscInt)(n*PetscRealPart(a)); 3251 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3252 } 3253 } 3254 } else { 3255 for (i=0; i<aij->nz; i++) {ierr = PetscRandomGetValue(rctx,aij->a+i);CHKERRQ(ierr);} 3256 } 3257 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3258 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3259 PetscFunctionReturn(0); 3260 } 3261 3262 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3263 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx) 3264 { 3265 PetscErrorCode ierr; 3266 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3267 PetscScalar a; 3268 PetscInt m,n,i,j,col,nskip; 3269 3270 PetscFunctionBegin; 3271 nskip = high - low; 3272 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3273 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3274 for (i=0; i<m; i++) { 3275 for (j=0; j<aij->imax[i]; j++) { 3276 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3277 col = (PetscInt)(n*PetscRealPart(a)); 3278 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3279 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3280 } 3281 } 3282 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3283 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3284 PetscFunctionReturn(0); 3285 } 3286 3287 3288 /* -------------------------------------------------------------------*/ 3289 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3290 MatGetRow_SeqAIJ, 3291 MatRestoreRow_SeqAIJ, 3292 MatMult_SeqAIJ, 3293 /* 4*/ MatMultAdd_SeqAIJ, 3294 MatMultTranspose_SeqAIJ, 3295 MatMultTransposeAdd_SeqAIJ, 3296 0, 3297 0, 3298 0, 3299 /* 10*/ 0, 3300 MatLUFactor_SeqAIJ, 3301 0, 3302 MatSOR_SeqAIJ, 3303 MatTranspose_SeqAIJ, 3304 /*1 5*/ MatGetInfo_SeqAIJ, 3305 MatEqual_SeqAIJ, 3306 MatGetDiagonal_SeqAIJ, 3307 MatDiagonalScale_SeqAIJ, 3308 MatNorm_SeqAIJ, 3309 /* 20*/ 0, 3310 MatAssemblyEnd_SeqAIJ, 3311 MatSetOption_SeqAIJ, 3312 MatZeroEntries_SeqAIJ, 3313 /* 24*/ MatZeroRows_SeqAIJ, 3314 0, 3315 0, 3316 0, 3317 0, 3318 /* 29*/ MatSetUp_SeqAIJ, 3319 0, 3320 0, 3321 0, 3322 0, 3323 /* 34*/ MatDuplicate_SeqAIJ, 3324 0, 3325 0, 3326 MatILUFactor_SeqAIJ, 3327 0, 3328 /* 39*/ MatAXPY_SeqAIJ, 3329 MatCreateSubMatrices_SeqAIJ, 3330 MatIncreaseOverlap_SeqAIJ, 3331 MatGetValues_SeqAIJ, 3332 MatCopy_SeqAIJ, 3333 /* 44*/ MatGetRowMax_SeqAIJ, 3334 MatScale_SeqAIJ, 3335 MatShift_SeqAIJ, 3336 MatDiagonalSet_SeqAIJ, 3337 MatZeroRowsColumns_SeqAIJ, 3338 /* 49*/ MatSetRandom_SeqAIJ, 3339 MatGetRowIJ_SeqAIJ, 3340 MatRestoreRowIJ_SeqAIJ, 3341 MatGetColumnIJ_SeqAIJ, 3342 MatRestoreColumnIJ_SeqAIJ, 3343 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3344 0, 3345 0, 3346 MatPermute_SeqAIJ, 3347 0, 3348 /* 59*/ 0, 3349 MatDestroy_SeqAIJ, 3350 MatView_SeqAIJ, 3351 0, 3352 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3353 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3354 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3355 0, 3356 0, 3357 0, 3358 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3359 MatGetRowMinAbs_SeqAIJ, 3360 0, 3361 0, 3362 0, 3363 /* 74*/ 0, 3364 MatFDColoringApply_AIJ, 3365 0, 3366 0, 3367 0, 3368 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3369 0, 3370 0, 3371 0, 3372 MatLoad_SeqAIJ, 3373 /* 84*/ MatIsSymmetric_SeqAIJ, 3374 MatIsHermitian_SeqAIJ, 3375 0, 3376 0, 3377 0, 3378 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3379 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3380 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3381 MatPtAP_SeqAIJ_SeqAIJ, 3382 MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy, 3383 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3384 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3385 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3386 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3387 0, 3388 /* 99*/ 0, 3389 0, 3390 0, 3391 MatConjugate_SeqAIJ, 3392 0, 3393 /*104*/ MatSetValuesRow_SeqAIJ, 3394 MatRealPart_SeqAIJ, 3395 MatImaginaryPart_SeqAIJ, 3396 0, 3397 0, 3398 /*109*/ MatMatSolve_SeqAIJ, 3399 0, 3400 MatGetRowMin_SeqAIJ, 3401 0, 3402 MatMissingDiagonal_SeqAIJ, 3403 /*114*/ 0, 3404 0, 3405 0, 3406 0, 3407 0, 3408 /*119*/ 0, 3409 0, 3410 0, 3411 0, 3412 MatGetMultiProcBlock_SeqAIJ, 3413 /*124*/ MatFindNonzeroRows_SeqAIJ, 3414 MatGetColumnNorms_SeqAIJ, 3415 MatInvertBlockDiagonal_SeqAIJ, 3416 MatInvertVariableBlockDiagonal_SeqAIJ, 3417 0, 3418 /*129*/ 0, 3419 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3420 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3421 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3422 MatTransposeColoringCreate_SeqAIJ, 3423 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3424 MatTransColoringApplyDenToSp_SeqAIJ, 3425 MatRARt_SeqAIJ_SeqAIJ, 3426 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3427 MatRARtNumeric_SeqAIJ_SeqAIJ, 3428 /*139*/0, 3429 0, 3430 0, 3431 MatFDColoringSetUp_SeqXAIJ, 3432 MatFindOffBlockDiagonalEntries_SeqAIJ, 3433 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3434 MatDestroySubMatrices_SeqAIJ 3435 }; 3436 3437 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3438 { 3439 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3440 PetscInt i,nz,n; 3441 3442 PetscFunctionBegin; 3443 nz = aij->maxnz; 3444 n = mat->rmap->n; 3445 for (i=0; i<nz; i++) { 3446 aij->j[i] = indices[i]; 3447 } 3448 aij->nz = nz; 3449 for (i=0; i<n; i++) { 3450 aij->ilen[i] = aij->imax[i]; 3451 } 3452 PetscFunctionReturn(0); 3453 } 3454 3455 /* 3456 * When a sparse matrix has many zero columns, we should compact them out to save the space 3457 * This happens in MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3458 * */ 3459 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3460 { 3461 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3462 PetscTable gid1_lid1; 3463 PetscTablePosition tpos; 3464 PetscInt gid,lid,i,j,ncols,ec; 3465 PetscInt *garray; 3466 PetscErrorCode ierr; 3467 3468 PetscFunctionBegin; 3469 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3470 PetscValidPointer(mapping,2); 3471 /* use a table */ 3472 ierr = PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1);CHKERRQ(ierr); 3473 ec = 0; 3474 for (i=0; i<mat->rmap->n; i++) { 3475 ncols = aij->i[i+1] - aij->i[i]; 3476 for (j=0; j<ncols; j++) { 3477 PetscInt data,gid1 = aij->j[aij->i[i] + j] + 1; 3478 ierr = PetscTableFind(gid1_lid1,gid1,&data);CHKERRQ(ierr); 3479 if (!data) { 3480 /* one based table */ 3481 ierr = PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);CHKERRQ(ierr); 3482 } 3483 } 3484 } 3485 /* form array of columns we need */ 3486 ierr = PetscMalloc1(ec+1,&garray);CHKERRQ(ierr); 3487 ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr); 3488 while (tpos) { 3489 ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);CHKERRQ(ierr); 3490 gid--; 3491 lid--; 3492 garray[lid] = gid; 3493 } 3494 ierr = PetscSortInt(ec,garray);CHKERRQ(ierr); /* sort, and rebuild */ 3495 ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr); 3496 for (i=0; i<ec; i++) { 3497 ierr = PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);CHKERRQ(ierr); 3498 } 3499 /* compact out the extra columns in B */ 3500 for (i=0; i<mat->rmap->n; i++) { 3501 ncols = aij->i[i+1] - aij->i[i]; 3502 for (j=0; j<ncols; j++) { 3503 PetscInt gid1 = aij->j[aij->i[i] + j] + 1; 3504 ierr = PetscTableFind(gid1_lid1,gid1,&lid);CHKERRQ(ierr); 3505 lid--; 3506 aij->j[aij->i[i] + j] = lid; 3507 } 3508 } 3509 mat->cmap->n = mat->cmap->N = ec; 3510 mat->cmap->bs = 1; 3511 3512 ierr = PetscTableDestroy(&gid1_lid1);CHKERRQ(ierr); 3513 ierr = PetscLayoutSetUp((mat->cmap));CHKERRQ(ierr); 3514 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping);CHKERRQ(ierr); 3515 ierr = ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH);CHKERRQ(ierr); 3516 PetscFunctionReturn(0); 3517 } 3518 3519 /*@ 3520 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3521 in the matrix. 3522 3523 Input Parameters: 3524 + mat - the SeqAIJ matrix 3525 - indices - the column indices 3526 3527 Level: advanced 3528 3529 Notes: 3530 This can be called if you have precomputed the nonzero structure of the 3531 matrix and want to provide it to the matrix object to improve the performance 3532 of the MatSetValues() operation. 3533 3534 You MUST have set the correct numbers of nonzeros per row in the call to 3535 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3536 3537 MUST be called before any calls to MatSetValues(); 3538 3539 The indices should start with zero, not one. 3540 3541 @*/ 3542 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3543 { 3544 PetscErrorCode ierr; 3545 3546 PetscFunctionBegin; 3547 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3548 PetscValidPointer(indices,2); 3549 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3550 PetscFunctionReturn(0); 3551 } 3552 3553 /* ----------------------------------------------------------------------------------------*/ 3554 3555 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3556 { 3557 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3558 PetscErrorCode ierr; 3559 size_t nz = aij->i[mat->rmap->n]; 3560 3561 PetscFunctionBegin; 3562 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3563 3564 /* allocate space for values if not already there */ 3565 if (!aij->saved_values) { 3566 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3567 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3568 } 3569 3570 /* copy values over */ 3571 ierr = PetscArraycpy(aij->saved_values,aij->a,nz);CHKERRQ(ierr); 3572 PetscFunctionReturn(0); 3573 } 3574 3575 /*@ 3576 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3577 example, reuse of the linear part of a Jacobian, while recomputing the 3578 nonlinear portion. 3579 3580 Collect on Mat 3581 3582 Input Parameters: 3583 . mat - the matrix (currently only AIJ matrices support this option) 3584 3585 Level: advanced 3586 3587 Common Usage, with SNESSolve(): 3588 $ Create Jacobian matrix 3589 $ Set linear terms into matrix 3590 $ Apply boundary conditions to matrix, at this time matrix must have 3591 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3592 $ boundary conditions again will not change the nonzero structure 3593 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3594 $ ierr = MatStoreValues(mat); 3595 $ Call SNESSetJacobian() with matrix 3596 $ In your Jacobian routine 3597 $ ierr = MatRetrieveValues(mat); 3598 $ Set nonlinear terms in matrix 3599 3600 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3601 $ // build linear portion of Jacobian 3602 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3603 $ ierr = MatStoreValues(mat); 3604 $ loop over nonlinear iterations 3605 $ ierr = MatRetrieveValues(mat); 3606 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3607 $ // call MatAssemblyBegin/End() on matrix 3608 $ Solve linear system with Jacobian 3609 $ endloop 3610 3611 Notes: 3612 Matrix must already be assemblied before calling this routine 3613 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3614 calling this routine. 3615 3616 When this is called multiple times it overwrites the previous set of stored values 3617 and does not allocated additional space. 3618 3619 .seealso: MatRetrieveValues() 3620 3621 @*/ 3622 PetscErrorCode MatStoreValues(Mat mat) 3623 { 3624 PetscErrorCode ierr; 3625 3626 PetscFunctionBegin; 3627 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3628 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3629 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3630 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3631 PetscFunctionReturn(0); 3632 } 3633 3634 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3635 { 3636 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3637 PetscErrorCode ierr; 3638 PetscInt nz = aij->i[mat->rmap->n]; 3639 3640 PetscFunctionBegin; 3641 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3642 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3643 /* copy values over */ 3644 ierr = PetscArraycpy(aij->a,aij->saved_values,nz);CHKERRQ(ierr); 3645 PetscFunctionReturn(0); 3646 } 3647 3648 /*@ 3649 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3650 example, reuse of the linear part of a Jacobian, while recomputing the 3651 nonlinear portion. 3652 3653 Collect on Mat 3654 3655 Input Parameters: 3656 . mat - the matrix (currently only AIJ matrices support this option) 3657 3658 Level: advanced 3659 3660 .seealso: MatStoreValues() 3661 3662 @*/ 3663 PetscErrorCode MatRetrieveValues(Mat mat) 3664 { 3665 PetscErrorCode ierr; 3666 3667 PetscFunctionBegin; 3668 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3669 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3670 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3671 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3672 PetscFunctionReturn(0); 3673 } 3674 3675 3676 /* --------------------------------------------------------------------------------*/ 3677 /*@C 3678 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3679 (the default parallel PETSc format). For good matrix assembly performance 3680 the user should preallocate the matrix storage by setting the parameter nz 3681 (or the array nnz). By setting these parameters accurately, performance 3682 during matrix assembly can be increased by more than a factor of 50. 3683 3684 Collective 3685 3686 Input Parameters: 3687 + comm - MPI communicator, set to PETSC_COMM_SELF 3688 . m - number of rows 3689 . n - number of columns 3690 . nz - number of nonzeros per row (same for all rows) 3691 - nnz - array containing the number of nonzeros in the various rows 3692 (possibly different for each row) or NULL 3693 3694 Output Parameter: 3695 . A - the matrix 3696 3697 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3698 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3699 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3700 3701 Notes: 3702 If nnz is given then nz is ignored 3703 3704 The AIJ format (also called the Yale sparse matrix format or 3705 compressed row storage), is fully compatible with standard Fortran 77 3706 storage. That is, the stored row and column indices can begin at 3707 either one (as in Fortran) or zero. See the users' manual for details. 3708 3709 Specify the preallocated storage with either nz or nnz (not both). 3710 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3711 allocation. For large problems you MUST preallocate memory or you 3712 will get TERRIBLE performance, see the users' manual chapter on matrices. 3713 3714 By default, this format uses inodes (identical nodes) when possible, to 3715 improve numerical efficiency of matrix-vector products and solves. We 3716 search for consecutive rows with the same nonzero structure, thereby 3717 reusing matrix information to achieve increased efficiency. 3718 3719 Options Database Keys: 3720 + -mat_no_inode - Do not use inodes 3721 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3722 3723 Level: intermediate 3724 3725 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3726 3727 @*/ 3728 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3729 { 3730 PetscErrorCode ierr; 3731 3732 PetscFunctionBegin; 3733 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3734 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3735 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3736 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3737 PetscFunctionReturn(0); 3738 } 3739 3740 /*@C 3741 MatSeqAIJSetPreallocation - For good matrix assembly performance 3742 the user should preallocate the matrix storage by setting the parameter nz 3743 (or the array nnz). By setting these parameters accurately, performance 3744 during matrix assembly can be increased by more than a factor of 50. 3745 3746 Collective 3747 3748 Input Parameters: 3749 + B - The matrix 3750 . nz - number of nonzeros per row (same for all rows) 3751 - nnz - array containing the number of nonzeros in the various rows 3752 (possibly different for each row) or NULL 3753 3754 Notes: 3755 If nnz is given then nz is ignored 3756 3757 The AIJ format (also called the Yale sparse matrix format or 3758 compressed row storage), is fully compatible with standard Fortran 77 3759 storage. That is, the stored row and column indices can begin at 3760 either one (as in Fortran) or zero. See the users' manual for details. 3761 3762 Specify the preallocated storage with either nz or nnz (not both). 3763 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3764 allocation. For large problems you MUST preallocate memory or you 3765 will get TERRIBLE performance, see the users' manual chapter on matrices. 3766 3767 You can call MatGetInfo() to get information on how effective the preallocation was; 3768 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3769 You can also run with the option -info and look for messages with the string 3770 malloc in them to see if additional memory allocation was needed. 3771 3772 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3773 entries or columns indices 3774 3775 By default, this format uses inodes (identical nodes) when possible, to 3776 improve numerical efficiency of matrix-vector products and solves. We 3777 search for consecutive rows with the same nonzero structure, thereby 3778 reusing matrix information to achieve increased efficiency. 3779 3780 Options Database Keys: 3781 + -mat_no_inode - Do not use inodes 3782 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3783 3784 Level: intermediate 3785 3786 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3787 3788 @*/ 3789 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3790 { 3791 PetscErrorCode ierr; 3792 3793 PetscFunctionBegin; 3794 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3795 PetscValidType(B,1); 3796 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3797 PetscFunctionReturn(0); 3798 } 3799 3800 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3801 { 3802 Mat_SeqAIJ *b; 3803 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3804 PetscErrorCode ierr; 3805 PetscInt i; 3806 3807 PetscFunctionBegin; 3808 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3809 if (nz == MAT_SKIP_ALLOCATION) { 3810 skipallocation = PETSC_TRUE; 3811 nz = 0; 3812 } 3813 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3814 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3815 3816 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3817 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3818 #if defined(PETSC_USE_DEBUG) 3819 if (nnz) { 3820 for (i=0; i<B->rmap->n; i++) { 3821 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]); 3822 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); 3823 } 3824 } 3825 #endif 3826 3827 B->preallocated = PETSC_TRUE; 3828 3829 b = (Mat_SeqAIJ*)B->data; 3830 3831 if (!skipallocation) { 3832 if (!b->imax) { 3833 ierr = PetscMalloc1(B->rmap->n,&b->imax);CHKERRQ(ierr); 3834 ierr = PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3835 } 3836 if (!b->ilen) { 3837 /* b->ilen will count nonzeros in each row so far. */ 3838 ierr = PetscCalloc1(B->rmap->n,&b->ilen);CHKERRQ(ierr); 3839 ierr = PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3840 } else { 3841 ierr = PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3842 } 3843 if (!b->ipre) { 3844 ierr = PetscMalloc1(B->rmap->n,&b->ipre);CHKERRQ(ierr); 3845 ierr = PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3846 } 3847 if (!nnz) { 3848 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3849 else if (nz < 0) nz = 1; 3850 nz = PetscMin(nz,B->cmap->n); 3851 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3852 nz = nz*B->rmap->n; 3853 } else { 3854 nz = 0; 3855 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3856 } 3857 3858 /* allocate the matrix space */ 3859 /* FIXME: should B's old memory be unlogged? */ 3860 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3861 if (B->structure_only) { 3862 ierr = PetscMalloc1(nz,&b->j);CHKERRQ(ierr); 3863 ierr = PetscMalloc1(B->rmap->n+1,&b->i);CHKERRQ(ierr); 3864 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));CHKERRQ(ierr); 3865 } else { 3866 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3867 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3868 } 3869 b->i[0] = 0; 3870 for (i=1; i<B->rmap->n+1; i++) { 3871 b->i[i] = b->i[i-1] + b->imax[i-1]; 3872 } 3873 if (B->structure_only) { 3874 b->singlemalloc = PETSC_FALSE; 3875 b->free_a = PETSC_FALSE; 3876 } else { 3877 b->singlemalloc = PETSC_TRUE; 3878 b->free_a = PETSC_TRUE; 3879 } 3880 b->free_ij = PETSC_TRUE; 3881 } else { 3882 b->free_a = PETSC_FALSE; 3883 b->free_ij = PETSC_FALSE; 3884 } 3885 3886 if (b->ipre && nnz != b->ipre && b->imax) { 3887 /* reserve user-requested sparsity */ 3888 ierr = PetscArraycpy(b->ipre,b->imax,B->rmap->n);CHKERRQ(ierr); 3889 } 3890 3891 3892 b->nz = 0; 3893 b->maxnz = nz; 3894 B->info.nz_unneeded = (double)b->maxnz; 3895 if (realalloc) { 3896 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3897 } 3898 B->was_assembled = PETSC_FALSE; 3899 B->assembled = PETSC_FALSE; 3900 PetscFunctionReturn(0); 3901 } 3902 3903 3904 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 3905 { 3906 Mat_SeqAIJ *a; 3907 PetscInt i; 3908 PetscErrorCode ierr; 3909 3910 PetscFunctionBegin; 3911 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3912 a = (Mat_SeqAIJ*)A->data; 3913 /* if no saved info, we error out */ 3914 if (!a->ipre) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"No saved preallocation info \n"); 3915 3916 if (!a->i || !a->j || !a->a || !a->imax || !a->ilen) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"Memory info is incomplete, and can not reset preallocation \n"); 3917 3918 ierr = PetscArraycpy(a->imax,a->ipre,A->rmap->n);CHKERRQ(ierr); 3919 ierr = PetscArrayzero(a->ilen,A->rmap->n);CHKERRQ(ierr); 3920 a->i[0] = 0; 3921 for (i=1; i<A->rmap->n+1; i++) { 3922 a->i[i] = a->i[i-1] + a->imax[i-1]; 3923 } 3924 A->preallocated = PETSC_TRUE; 3925 a->nz = 0; 3926 a->maxnz = a->i[A->rmap->n]; 3927 A->info.nz_unneeded = (double)a->maxnz; 3928 A->was_assembled = PETSC_FALSE; 3929 A->assembled = PETSC_FALSE; 3930 PetscFunctionReturn(0); 3931 } 3932 3933 /*@ 3934 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3935 3936 Input Parameters: 3937 + B - the matrix 3938 . i - the indices into j for the start of each row (starts with zero) 3939 . j - the column indices for each row (starts with zero) these must be sorted for each row 3940 - v - optional values in the matrix 3941 3942 Level: developer 3943 3944 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3945 3946 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ 3947 @*/ 3948 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3949 { 3950 PetscErrorCode ierr; 3951 3952 PetscFunctionBegin; 3953 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3954 PetscValidType(B,1); 3955 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3956 PetscFunctionReturn(0); 3957 } 3958 3959 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3960 { 3961 PetscInt i; 3962 PetscInt m,n; 3963 PetscInt nz; 3964 PetscInt *nnz, nz_max = 0; 3965 PetscErrorCode ierr; 3966 3967 PetscFunctionBegin; 3968 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3969 3970 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3971 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3972 3973 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3974 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3975 for (i = 0; i < m; i++) { 3976 nz = Ii[i+1]- Ii[i]; 3977 nz_max = PetscMax(nz_max, nz); 3978 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3979 nnz[i] = nz; 3980 } 3981 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3982 ierr = PetscFree(nnz);CHKERRQ(ierr); 3983 3984 for (i = 0; i < m; i++) { 3985 ierr = MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES);CHKERRQ(ierr); 3986 } 3987 3988 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3989 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3990 3991 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3992 PetscFunctionReturn(0); 3993 } 3994 3995 #include <../src/mat/impls/dense/seq/dense.h> 3996 #include <petsc/private/kernels/petscaxpy.h> 3997 3998 /* 3999 Computes (B'*A')' since computing B*A directly is untenable 4000 4001 n p p 4002 ( ) ( ) ( ) 4003 m ( A ) * n ( B ) = m ( C ) 4004 ( ) ( ) ( ) 4005 4006 */ 4007 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 4008 { 4009 PetscErrorCode ierr; 4010 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 4011 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 4012 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 4013 PetscInt i,n,m,q,p; 4014 const PetscInt *ii,*idx; 4015 const PetscScalar *b,*a,*a_q; 4016 PetscScalar *c,*c_q; 4017 4018 PetscFunctionBegin; 4019 m = A->rmap->n; 4020 n = A->cmap->n; 4021 p = B->cmap->n; 4022 a = sub_a->v; 4023 b = sub_b->a; 4024 c = sub_c->v; 4025 ierr = PetscArrayzero(c,m*p);CHKERRQ(ierr); 4026 4027 ii = sub_b->i; 4028 idx = sub_b->j; 4029 for (i=0; i<n; i++) { 4030 q = ii[i+1] - ii[i]; 4031 while (q-->0) { 4032 c_q = c + m*(*idx); 4033 a_q = a + m*i; 4034 PetscKernelAXPY(c_q,*b,a_q,m); 4035 idx++; 4036 b++; 4037 } 4038 } 4039 PetscFunctionReturn(0); 4040 } 4041 4042 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 4043 { 4044 PetscErrorCode ierr; 4045 PetscInt m=A->rmap->n,n=B->cmap->n; 4046 Mat Cmat; 4047 4048 PetscFunctionBegin; 4049 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); 4050 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 4051 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 4052 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 4053 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 4054 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 4055 4056 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4057 4058 *C = Cmat; 4059 PetscFunctionReturn(0); 4060 } 4061 4062 /* ----------------------------------------------------------------*/ 4063 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 4064 { 4065 PetscErrorCode ierr; 4066 4067 PetscFunctionBegin; 4068 if (scall == MAT_INITIAL_MATRIX) { 4069 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 4070 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 4071 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 4072 } 4073 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 4074 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 4075 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 4076 PetscFunctionReturn(0); 4077 } 4078 4079 4080 /*MC 4081 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4082 based on compressed sparse row format. 4083 4084 Options Database Keys: 4085 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4086 4087 Level: beginner 4088 4089 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 4090 M*/ 4091 4092 /*MC 4093 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4094 4095 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4096 and MATMPIAIJ otherwise. As a result, for single process communicators, 4097 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 4098 for communicators controlling multiple processes. It is recommended that you call both of 4099 the above preallocation routines for simplicity. 4100 4101 Options Database Keys: 4102 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4103 4104 Developer Notes: 4105 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4106 enough exist. 4107 4108 Level: beginner 4109 4110 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 4111 M*/ 4112 4113 /*MC 4114 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4115 4116 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4117 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4118 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4119 for communicators controlling multiple processes. It is recommended that you call both of 4120 the above preallocation routines for simplicity. 4121 4122 Options Database Keys: 4123 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4124 4125 Level: beginner 4126 4127 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 4128 M*/ 4129 4130 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4131 #if defined(PETSC_HAVE_ELEMENTAL) 4132 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4133 #endif 4134 #if defined(PETSC_HAVE_HYPRE) 4135 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4136 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 4137 #endif 4138 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 4139 4140 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4141 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4142 PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*); 4143 4144 /*@C 4145 MatSeqAIJGetArray - gives access to the array where the data for a MATSEQAIJ matrix is stored 4146 4147 Not Collective 4148 4149 Input Parameter: 4150 . mat - a MATSEQAIJ matrix 4151 4152 Output Parameter: 4153 . array - pointer to the data 4154 4155 Level: intermediate 4156 4157 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4158 @*/ 4159 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4160 { 4161 PetscErrorCode ierr; 4162 4163 PetscFunctionBegin; 4164 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4165 PetscFunctionReturn(0); 4166 } 4167 4168 /*@C 4169 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4170 4171 Not Collective 4172 4173 Input Parameter: 4174 . mat - a MATSEQAIJ matrix 4175 4176 Output Parameter: 4177 . nz - the maximum number of nonzeros in any row 4178 4179 Level: intermediate 4180 4181 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4182 @*/ 4183 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4184 { 4185 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4186 4187 PetscFunctionBegin; 4188 *nz = aij->rmax; 4189 PetscFunctionReturn(0); 4190 } 4191 4192 /*@C 4193 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4194 4195 Not Collective 4196 4197 Input Parameters: 4198 . mat - a MATSEQAIJ matrix 4199 . array - pointer to the data 4200 4201 Level: intermediate 4202 4203 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4204 @*/ 4205 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4206 { 4207 PetscErrorCode ierr; 4208 4209 PetscFunctionBegin; 4210 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 4211 PetscFunctionReturn(0); 4212 } 4213 4214 #if defined(PETSC_HAVE_CUDA) 4215 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat); 4216 #endif 4217 4218 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4219 { 4220 Mat_SeqAIJ *b; 4221 PetscErrorCode ierr; 4222 PetscMPIInt size; 4223 4224 PetscFunctionBegin; 4225 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4226 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4227 4228 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4229 4230 B->data = (void*)b; 4231 4232 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4233 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4234 4235 b->row = 0; 4236 b->col = 0; 4237 b->icol = 0; 4238 b->reallocs = 0; 4239 b->ignorezeroentries = PETSC_FALSE; 4240 b->roworiented = PETSC_TRUE; 4241 b->nonew = 0; 4242 b->diag = 0; 4243 b->solve_work = 0; 4244 B->spptr = 0; 4245 b->saved_values = 0; 4246 b->idiag = 0; 4247 b->mdiag = 0; 4248 b->ssor_work = 0; 4249 b->omega = 1.0; 4250 b->fshift = 0.0; 4251 b->idiagvalid = PETSC_FALSE; 4252 b->ibdiagvalid = PETSC_FALSE; 4253 b->keepnonzeropattern = PETSC_FALSE; 4254 4255 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4256 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4257 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4258 4259 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4260 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4261 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4262 #endif 4263 4264 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4265 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4266 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4267 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4268 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4269 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4270 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL);CHKERRQ(ierr); 4271 #if defined(PETSC_HAVE_MKL_SPARSE) 4272 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4273 #endif 4274 #if defined(PETSC_HAVE_CUDA) 4275 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE);CHKERRQ(ierr); 4276 #endif 4277 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4278 #if defined(PETSC_HAVE_ELEMENTAL) 4279 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4280 #endif 4281 #if defined(PETSC_HAVE_HYPRE) 4282 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 4283 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 4284 #endif 4285 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 4286 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);CHKERRQ(ierr); 4287 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS);CHKERRQ(ierr); 4288 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4289 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4290 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4291 ierr = PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);CHKERRQ(ierr); 4292 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4293 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4294 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4295 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4296 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4297 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_seqaij_C",MatPtAP_IS_XAIJ);CHKERRQ(ierr); 4298 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4299 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4300 ierr = MatSeqAIJSetTypeFromOptions(B);CHKERRQ(ierr); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4301 PetscFunctionReturn(0); 4302 } 4303 4304 /* 4305 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4306 */ 4307 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4308 { 4309 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4310 PetscErrorCode ierr; 4311 PetscInt m = A->rmap->n,i; 4312 4313 PetscFunctionBegin; 4314 c = (Mat_SeqAIJ*)C->data; 4315 4316 C->factortype = A->factortype; 4317 c->row = 0; 4318 c->col = 0; 4319 c->icol = 0; 4320 c->reallocs = 0; 4321 4322 C->assembled = PETSC_TRUE; 4323 4324 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4325 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4326 4327 ierr = PetscMalloc1(m,&c->imax);CHKERRQ(ierr); 4328 ierr = PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt));CHKERRQ(ierr); 4329 ierr = PetscMalloc1(m,&c->ilen);CHKERRQ(ierr); 4330 ierr = PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt));CHKERRQ(ierr); 4331 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4332 4333 /* allocate the matrix space */ 4334 if (mallocmatspace) { 4335 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4336 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4337 4338 c->singlemalloc = PETSC_TRUE; 4339 4340 ierr = PetscArraycpy(c->i,a->i,m+1);CHKERRQ(ierr); 4341 if (m > 0) { 4342 ierr = PetscArraycpy(c->j,a->j,a->i[m]);CHKERRQ(ierr); 4343 if (cpvalues == MAT_COPY_VALUES) { 4344 ierr = PetscArraycpy(c->a,a->a,a->i[m]);CHKERRQ(ierr); 4345 } else { 4346 ierr = PetscArrayzero(c->a,a->i[m]);CHKERRQ(ierr); 4347 } 4348 } 4349 } 4350 4351 c->ignorezeroentries = a->ignorezeroentries; 4352 c->roworiented = a->roworiented; 4353 c->nonew = a->nonew; 4354 if (a->diag) { 4355 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4356 ierr = PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt));CHKERRQ(ierr); 4357 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4358 } else c->diag = NULL; 4359 4360 c->solve_work = 0; 4361 c->saved_values = 0; 4362 c->idiag = 0; 4363 c->ssor_work = 0; 4364 c->keepnonzeropattern = a->keepnonzeropattern; 4365 c->free_a = PETSC_TRUE; 4366 c->free_ij = PETSC_TRUE; 4367 4368 c->rmax = a->rmax; 4369 c->nz = a->nz; 4370 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4371 C->preallocated = PETSC_TRUE; 4372 4373 c->compressedrow.use = a->compressedrow.use; 4374 c->compressedrow.nrows = a->compressedrow.nrows; 4375 if (a->compressedrow.use) { 4376 i = a->compressedrow.nrows; 4377 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4378 ierr = PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1);CHKERRQ(ierr); 4379 ierr = PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i);CHKERRQ(ierr); 4380 } else { 4381 c->compressedrow.use = PETSC_FALSE; 4382 c->compressedrow.i = NULL; 4383 c->compressedrow.rindex = NULL; 4384 } 4385 c->nonzerorowcnt = a->nonzerorowcnt; 4386 C->nonzerostate = A->nonzerostate; 4387 4388 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4389 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4390 PetscFunctionReturn(0); 4391 } 4392 4393 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4394 { 4395 PetscErrorCode ierr; 4396 4397 PetscFunctionBegin; 4398 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4399 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4400 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4401 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4402 } 4403 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4404 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4405 PetscFunctionReturn(0); 4406 } 4407 4408 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4409 { 4410 PetscBool isbinary, ishdf5; 4411 PetscErrorCode ierr; 4412 4413 PetscFunctionBegin; 4414 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4415 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4416 /* force binary viewer to load .info file if it has not yet done so */ 4417 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4418 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 4419 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 4420 if (isbinary) { 4421 ierr = MatLoad_SeqAIJ_Binary(newMat,viewer);CHKERRQ(ierr); 4422 } else if (ishdf5) { 4423 #if defined(PETSC_HAVE_HDF5) 4424 ierr = MatLoad_AIJ_HDF5(newMat,viewer);CHKERRQ(ierr); 4425 #else 4426 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4427 #endif 4428 } else { 4429 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); 4430 } 4431 PetscFunctionReturn(0); 4432 } 4433 4434 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat newMat, PetscViewer viewer) 4435 { 4436 Mat_SeqAIJ *a; 4437 PetscErrorCode ierr; 4438 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4439 int fd; 4440 PetscMPIInt size; 4441 MPI_Comm comm; 4442 PetscInt bs = newMat->rmap->bs; 4443 4444 PetscFunctionBegin; 4445 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4446 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4447 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4448 4449 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4450 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4451 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4452 if (bs < 0) bs = 1; 4453 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4454 4455 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4456 ierr = PetscBinaryRead(fd,header,4,NULL,PETSC_INT);CHKERRQ(ierr); 4457 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4458 M = header[1]; N = header[2]; nz = header[3]; 4459 4460 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4461 4462 /* read in row lengths */ 4463 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4464 ierr = PetscBinaryRead(fd,rowlengths,M,NULL,PETSC_INT);CHKERRQ(ierr); 4465 4466 /* check if sum of rowlengths is same as nz */ 4467 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4468 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); 4469 4470 /* set global size if not set already*/ 4471 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4472 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4473 } else { 4474 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4475 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4476 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4477 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4478 } 4479 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); 4480 } 4481 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4482 a = (Mat_SeqAIJ*)newMat->data; 4483 4484 ierr = PetscBinaryRead(fd,a->j,nz,NULL,PETSC_INT);CHKERRQ(ierr); 4485 4486 /* read in nonzero values */ 4487 ierr = PetscBinaryRead(fd,a->a,nz,NULL,PETSC_SCALAR);CHKERRQ(ierr); 4488 4489 /* set matrix "i" values */ 4490 a->i[0] = 0; 4491 for (i=1; i<= M; i++) { 4492 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4493 a->ilen[i-1] = rowlengths[i-1]; 4494 } 4495 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4496 4497 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4498 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4499 PetscFunctionReturn(0); 4500 } 4501 4502 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4503 { 4504 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4505 PetscErrorCode ierr; 4506 #if defined(PETSC_USE_COMPLEX) 4507 PetscInt k; 4508 #endif 4509 4510 PetscFunctionBegin; 4511 /* If the matrix dimensions are not equal,or no of nonzeros */ 4512 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4513 *flg = PETSC_FALSE; 4514 PetscFunctionReturn(0); 4515 } 4516 4517 /* if the a->i are the same */ 4518 ierr = PetscArraycmp(a->i,b->i,A->rmap->n+1,flg);CHKERRQ(ierr); 4519 if (!*flg) PetscFunctionReturn(0); 4520 4521 /* if a->j are the same */ 4522 ierr = PetscArraycmp(a->j,b->j,a->nz,flg);CHKERRQ(ierr); 4523 if (!*flg) PetscFunctionReturn(0); 4524 4525 /* if a->a are the same */ 4526 #if defined(PETSC_USE_COMPLEX) 4527 for (k=0; k<a->nz; k++) { 4528 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4529 *flg = PETSC_FALSE; 4530 PetscFunctionReturn(0); 4531 } 4532 } 4533 #else 4534 ierr = PetscArraycmp(a->a,b->a,a->nz,flg);CHKERRQ(ierr); 4535 #endif 4536 PetscFunctionReturn(0); 4537 } 4538 4539 /*@ 4540 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4541 provided by the user. 4542 4543 Collective 4544 4545 Input Parameters: 4546 + comm - must be an MPI communicator of size 1 4547 . m - number of rows 4548 . n - number of columns 4549 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 4550 . j - column indices 4551 - a - matrix values 4552 4553 Output Parameter: 4554 . mat - the matrix 4555 4556 Level: intermediate 4557 4558 Notes: 4559 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4560 once the matrix is destroyed and not before 4561 4562 You cannot set new nonzero locations into this matrix, that will generate an error. 4563 4564 The i and j indices are 0 based 4565 4566 The format which is used for the sparse matrix input, is equivalent to a 4567 row-major ordering.. i.e for the following matrix, the input data expected is 4568 as shown 4569 4570 $ 1 0 0 4571 $ 2 0 3 4572 $ 4 5 6 4573 $ 4574 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4575 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4576 $ v = {1,2,3,4,5,6} [size = 6] 4577 4578 4579 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4580 4581 @*/ 4582 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4583 { 4584 PetscErrorCode ierr; 4585 PetscInt ii; 4586 Mat_SeqAIJ *aij; 4587 #if defined(PETSC_USE_DEBUG) 4588 PetscInt jj; 4589 #endif 4590 4591 PetscFunctionBegin; 4592 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4593 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4594 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4595 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4596 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4597 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4598 aij = (Mat_SeqAIJ*)(*mat)->data; 4599 ierr = PetscMalloc1(m,&aij->imax);CHKERRQ(ierr); 4600 ierr = PetscMalloc1(m,&aij->ilen);CHKERRQ(ierr); 4601 4602 aij->i = i; 4603 aij->j = j; 4604 aij->a = a; 4605 aij->singlemalloc = PETSC_FALSE; 4606 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4607 aij->free_a = PETSC_FALSE; 4608 aij->free_ij = PETSC_FALSE; 4609 4610 for (ii=0; ii<m; ii++) { 4611 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4612 #if defined(PETSC_USE_DEBUG) 4613 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]); 4614 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4615 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); 4616 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); 4617 } 4618 #endif 4619 } 4620 #if defined(PETSC_USE_DEBUG) 4621 for (ii=0; ii<aij->i[m]; ii++) { 4622 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4623 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]); 4624 } 4625 #endif 4626 4627 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4628 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4629 PetscFunctionReturn(0); 4630 } 4631 /*@C 4632 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4633 provided by the user. 4634 4635 Collective 4636 4637 Input Parameters: 4638 + comm - must be an MPI communicator of size 1 4639 . m - number of rows 4640 . n - number of columns 4641 . i - row indices 4642 . j - column indices 4643 . a - matrix values 4644 . nz - number of nonzeros 4645 - idx - 0 or 1 based 4646 4647 Output Parameter: 4648 . mat - the matrix 4649 4650 Level: intermediate 4651 4652 Notes: 4653 The i and j indices are 0 based 4654 4655 The format which is used for the sparse matrix input, is equivalent to a 4656 row-major ordering.. i.e for the following matrix, the input data expected is 4657 as shown: 4658 4659 1 0 0 4660 2 0 3 4661 4 5 6 4662 4663 i = {0,1,1,2,2,2} 4664 j = {0,0,2,0,1,2} 4665 v = {1,2,3,4,5,6} 4666 4667 4668 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4669 4670 @*/ 4671 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4672 { 4673 PetscErrorCode ierr; 4674 PetscInt ii, *nnz, one = 1,row,col; 4675 4676 4677 PetscFunctionBegin; 4678 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4679 for (ii = 0; ii < nz; ii++) { 4680 nnz[i[ii] - !!idx] += 1; 4681 } 4682 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4683 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4684 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4685 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4686 for (ii = 0; ii < nz; ii++) { 4687 if (idx) { 4688 row = i[ii] - 1; 4689 col = j[ii] - 1; 4690 } else { 4691 row = i[ii]; 4692 col = j[ii]; 4693 } 4694 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4695 } 4696 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4697 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4698 ierr = PetscFree(nnz);CHKERRQ(ierr); 4699 PetscFunctionReturn(0); 4700 } 4701 4702 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4703 { 4704 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4705 PetscErrorCode ierr; 4706 4707 PetscFunctionBegin; 4708 a->idiagvalid = PETSC_FALSE; 4709 a->ibdiagvalid = PETSC_FALSE; 4710 4711 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4712 PetscFunctionReturn(0); 4713 } 4714 4715 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4716 { 4717 PetscErrorCode ierr; 4718 PetscMPIInt size; 4719 4720 PetscFunctionBegin; 4721 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4722 if (size == 1) { 4723 if (scall == MAT_INITIAL_MATRIX) { 4724 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 4725 } else { 4726 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4727 } 4728 } else { 4729 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4730 } 4731 PetscFunctionReturn(0); 4732 } 4733 4734 /* 4735 Permute A into C's *local* index space using rowemb,colemb. 4736 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4737 of [0,m), colemb is in [0,n). 4738 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4739 */ 4740 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4741 { 4742 /* If making this function public, change the error returned in this function away from _PLIB. */ 4743 PetscErrorCode ierr; 4744 Mat_SeqAIJ *Baij; 4745 PetscBool seqaij; 4746 PetscInt m,n,*nz,i,j,count; 4747 PetscScalar v; 4748 const PetscInt *rowindices,*colindices; 4749 4750 PetscFunctionBegin; 4751 if (!B) PetscFunctionReturn(0); 4752 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4753 ierr = PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4754 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4755 if (rowemb) { 4756 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4757 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); 4758 } else { 4759 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4760 } 4761 if (colemb) { 4762 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4763 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); 4764 } else { 4765 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4766 } 4767 4768 Baij = (Mat_SeqAIJ*)(B->data); 4769 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4770 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4771 for (i=0; i<B->rmap->n; i++) { 4772 nz[i] = Baij->i[i+1] - Baij->i[i]; 4773 } 4774 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4775 ierr = PetscFree(nz);CHKERRQ(ierr); 4776 } 4777 if (pattern == SUBSET_NONZERO_PATTERN) { 4778 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4779 } 4780 count = 0; 4781 rowindices = NULL; 4782 colindices = NULL; 4783 if (rowemb) { 4784 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4785 } 4786 if (colemb) { 4787 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4788 } 4789 for (i=0; i<B->rmap->n; i++) { 4790 PetscInt row; 4791 row = i; 4792 if (rowindices) row = rowindices[i]; 4793 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4794 PetscInt col; 4795 col = Baij->j[count]; 4796 if (colindices) col = colindices[col]; 4797 v = Baij->a[count]; 4798 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4799 ++count; 4800 } 4801 } 4802 /* FIXME: set C's nonzerostate correctly. */ 4803 /* Assembly for C is necessary. */ 4804 C->preallocated = PETSC_TRUE; 4805 C->assembled = PETSC_TRUE; 4806 C->was_assembled = PETSC_FALSE; 4807 PetscFunctionReturn(0); 4808 } 4809 4810 PetscFunctionList MatSeqAIJList = NULL; 4811 4812 /*@C 4813 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 4814 4815 Collective on Mat 4816 4817 Input Parameters: 4818 + mat - the matrix object 4819 - matype - matrix type 4820 4821 Options Database Key: 4822 . -mat_seqai_type <method> - for example seqaijcrl 4823 4824 4825 Level: intermediate 4826 4827 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 4828 @*/ 4829 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 4830 { 4831 PetscErrorCode ierr,(*r)(Mat,MatType,MatReuse,Mat*); 4832 PetscBool sametype; 4833 4834 PetscFunctionBegin; 4835 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 4836 ierr = PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);CHKERRQ(ierr); 4837 if (sametype) PetscFunctionReturn(0); 4838 4839 ierr = PetscFunctionListFind(MatSeqAIJList,matype,&r);CHKERRQ(ierr); 4840 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 4841 ierr = (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);CHKERRQ(ierr); 4842 PetscFunctionReturn(0); 4843 } 4844 4845 4846 /*@C 4847 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 4848 4849 Not Collective 4850 4851 Input Parameters: 4852 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 4853 - function - routine to convert to subtype 4854 4855 Notes: 4856 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 4857 4858 4859 Then, your matrix can be chosen with the procedural interface at runtime via the option 4860 $ -mat_seqaij_type my_mat 4861 4862 Level: advanced 4863 4864 .seealso: MatSeqAIJRegisterAll() 4865 4866 4867 Level: advanced 4868 @*/ 4869 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 4870 { 4871 PetscErrorCode ierr; 4872 4873 PetscFunctionBegin; 4874 ierr = MatInitializePackage();CHKERRQ(ierr); 4875 ierr = PetscFunctionListAdd(&MatSeqAIJList,sname,function);CHKERRQ(ierr); 4876 PetscFunctionReturn(0); 4877 } 4878 4879 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 4880 4881 /*@C 4882 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 4883 4884 Not Collective 4885 4886 Level: advanced 4887 4888 Developers Note: CUSP and CUSPARSE do not yet support the MatConvert_SeqAIJ..() paradigm and thus cannot be registered here 4889 4890 .seealso: MatRegisterAll(), MatSeqAIJRegister() 4891 @*/ 4892 PetscErrorCode MatSeqAIJRegisterAll(void) 4893 { 4894 PetscErrorCode ierr; 4895 4896 PetscFunctionBegin; 4897 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 4898 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 4899 4900 ierr = MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4901 ierr = MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4902 ierr = MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL);CHKERRQ(ierr); 4903 #if defined(PETSC_HAVE_MKL_SPARSE) 4904 ierr = MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4905 #endif 4906 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 4907 ierr = MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr); 4908 #endif 4909 PetscFunctionReturn(0); 4910 } 4911 4912 /* 4913 Special version for direct calls from Fortran 4914 */ 4915 #include <petsc/private/fortranimpl.h> 4916 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4917 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4918 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4919 #define matsetvaluesseqaij_ matsetvaluesseqaij 4920 #endif 4921 4922 /* Change these macros so can be used in void function */ 4923 #undef CHKERRQ 4924 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4925 #undef SETERRQ2 4926 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4927 #undef SETERRQ3 4928 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4929 4930 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) 4931 { 4932 Mat A = *AA; 4933 PetscInt m = *mm, n = *nn; 4934 InsertMode is = *isis; 4935 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4936 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4937 PetscInt *imax,*ai,*ailen; 4938 PetscErrorCode ierr; 4939 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4940 MatScalar *ap,value,*aa; 4941 PetscBool ignorezeroentries = a->ignorezeroentries; 4942 PetscBool roworiented = a->roworiented; 4943 4944 PetscFunctionBegin; 4945 MatCheckPreallocated(A,1); 4946 imax = a->imax; 4947 ai = a->i; 4948 ailen = a->ilen; 4949 aj = a->j; 4950 aa = a->a; 4951 4952 for (k=0; k<m; k++) { /* loop over added rows */ 4953 row = im[k]; 4954 if (row < 0) continue; 4955 #if defined(PETSC_USE_DEBUG) 4956 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4957 #endif 4958 rp = aj + ai[row]; ap = aa + ai[row]; 4959 rmax = imax[row]; nrow = ailen[row]; 4960 low = 0; 4961 high = nrow; 4962 for (l=0; l<n; l++) { /* loop over added columns */ 4963 if (in[l] < 0) continue; 4964 #if defined(PETSC_USE_DEBUG) 4965 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4966 #endif 4967 col = in[l]; 4968 if (roworiented) value = v[l + k*n]; 4969 else value = v[k + l*m]; 4970 4971 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4972 4973 if (col <= lastcol) low = 0; 4974 else high = nrow; 4975 lastcol = col; 4976 while (high-low > 5) { 4977 t = (low+high)/2; 4978 if (rp[t] > col) high = t; 4979 else low = t; 4980 } 4981 for (i=low; i<high; i++) { 4982 if (rp[i] > col) break; 4983 if (rp[i] == col) { 4984 if (is == ADD_VALUES) ap[i] += value; 4985 else ap[i] = value; 4986 goto noinsert; 4987 } 4988 } 4989 if (value == 0.0 && ignorezeroentries) goto noinsert; 4990 if (nonew == 1) goto noinsert; 4991 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4992 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4993 N = nrow++ - 1; a->nz++; high++; 4994 /* shift up all the later entries in this row */ 4995 for (ii=N; ii>=i; ii--) { 4996 rp[ii+1] = rp[ii]; 4997 ap[ii+1] = ap[ii]; 4998 } 4999 rp[i] = col; 5000 ap[i] = value; 5001 A->nonzerostate++; 5002 noinsert:; 5003 low = i + 1; 5004 } 5005 ailen[row] = nrow; 5006 } 5007 PetscFunctionReturnVoid(); 5008 } 5009