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