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