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