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,d = 0; 1741 PetscErrorCode ierr; 1742 const PetscScalar *xx; 1743 PetscScalar *bb; 1744 PetscBool missing; 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 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1759 if (a->keepnonzeropattern) { 1760 for (i=0; i<N; i++) { 1761 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1762 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1763 } 1764 if (diag != 0.0) { 1765 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1766 for (i=0; i<N; i++) { 1767 a->a[a->diag[rows[i]]] = diag; 1768 } 1769 } 1770 } else { 1771 if (diag != 0.0) { 1772 for (i=0; i<N; i++) { 1773 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1774 if (a->ilen[rows[i]] > 0) { 1775 a->ilen[rows[i]] = 1; 1776 a->a[a->i[rows[i]]] = diag; 1777 a->j[a->i[rows[i]]] = rows[i]; 1778 } else { /* in case row was completely empty */ 1779 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1780 } 1781 } 1782 } else { 1783 for (i=0; i<N; i++) { 1784 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1785 a->ilen[rows[i]] = 0; 1786 } 1787 } 1788 A->nonzerostate++; 1789 } 1790 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1791 PetscFunctionReturn(0); 1792 } 1793 1794 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1795 { 1796 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1797 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1798 PetscErrorCode ierr; 1799 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 1800 const PetscScalar *xx; 1801 PetscScalar *bb; 1802 1803 PetscFunctionBegin; 1804 if (x && b) { 1805 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1806 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1807 vecs = PETSC_TRUE; 1808 } 1809 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1810 for (i=0; i<N; i++) { 1811 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1812 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1813 1814 zeroed[rows[i]] = PETSC_TRUE; 1815 } 1816 for (i=0; i<A->rmap->n; i++) { 1817 if (!zeroed[i]) { 1818 for (j=a->i[i]; j<a->i[i+1]; j++) { 1819 if (zeroed[a->j[j]]) { 1820 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 1821 a->a[j] = 0.0; 1822 } 1823 } 1824 } else if (vecs) bb[i] = diag*xx[i]; 1825 } 1826 if (x && b) { 1827 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1828 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1829 } 1830 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1831 if (diag != 0.0) { 1832 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1833 if (missing) { 1834 if (a->nonew) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1835 else { 1836 for (i=0; i<N; i++) { 1837 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1838 } 1839 } 1840 } else { 1841 for (i=0; i<N; i++) { 1842 a->a[a->diag[rows[i]]] = diag; 1843 } 1844 } 1845 } 1846 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1847 PetscFunctionReturn(0); 1848 } 1849 1850 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1851 { 1852 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1853 PetscInt *itmp; 1854 1855 PetscFunctionBegin; 1856 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1857 1858 *nz = a->i[row+1] - a->i[row]; 1859 if (v) *v = a->a + a->i[row]; 1860 if (idx) { 1861 itmp = a->j + a->i[row]; 1862 if (*nz) *idx = itmp; 1863 else *idx = 0; 1864 } 1865 PetscFunctionReturn(0); 1866 } 1867 1868 /* remove this function? */ 1869 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1870 { 1871 PetscFunctionBegin; 1872 PetscFunctionReturn(0); 1873 } 1874 1875 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1876 { 1877 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1878 MatScalar *v = a->a; 1879 PetscReal sum = 0.0; 1880 PetscErrorCode ierr; 1881 PetscInt i,j; 1882 1883 PetscFunctionBegin; 1884 if (type == NORM_FROBENIUS) { 1885 for (i=0; i<a->nz; i++) { 1886 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1887 } 1888 *nrm = PetscSqrtReal(sum); 1889 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 1890 } else if (type == NORM_1) { 1891 PetscReal *tmp; 1892 PetscInt *jj = a->j; 1893 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 1894 *nrm = 0.0; 1895 for (j=0; j<a->nz; j++) { 1896 tmp[*jj++] += PetscAbsScalar(*v); v++; 1897 } 1898 for (j=0; j<A->cmap->n; j++) { 1899 if (tmp[j] > *nrm) *nrm = tmp[j]; 1900 } 1901 ierr = PetscFree(tmp);CHKERRQ(ierr); 1902 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1903 } else if (type == NORM_INFINITY) { 1904 *nrm = 0.0; 1905 for (j=0; j<A->rmap->n; j++) { 1906 v = a->a + a->i[j]; 1907 sum = 0.0; 1908 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1909 sum += PetscAbsScalar(*v); v++; 1910 } 1911 if (sum > *nrm) *nrm = sum; 1912 } 1913 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1914 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1915 PetscFunctionReturn(0); 1916 } 1917 1918 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 1919 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 1920 { 1921 PetscErrorCode ierr; 1922 PetscInt i,j,anzj; 1923 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 1924 PetscInt an=A->cmap->N,am=A->rmap->N; 1925 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 1926 1927 PetscFunctionBegin; 1928 /* Allocate space for symbolic transpose info and work array */ 1929 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 1930 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 1931 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 1932 1933 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 1934 /* Note: offset by 1 for fast conversion into csr format. */ 1935 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 1936 /* Form ati for csr format of A^T. */ 1937 for (i=0;i<an;i++) ati[i+1] += ati[i]; 1938 1939 /* Copy ati into atfill so we have locations of the next free space in atj */ 1940 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 1941 1942 /* Walk through A row-wise and mark nonzero entries of A^T. */ 1943 for (i=0;i<am;i++) { 1944 anzj = ai[i+1] - ai[i]; 1945 for (j=0;j<anzj;j++) { 1946 atj[atfill[*aj]] = i; 1947 atfill[*aj++] += 1; 1948 } 1949 } 1950 1951 /* Clean up temporary space and complete requests. */ 1952 ierr = PetscFree(atfill);CHKERRQ(ierr); 1953 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 1954 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 1955 1956 b = (Mat_SeqAIJ*)((*B)->data); 1957 b->free_a = PETSC_FALSE; 1958 b->free_ij = PETSC_TRUE; 1959 b->nonew = 0; 1960 PetscFunctionReturn(0); 1961 } 1962 1963 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1964 { 1965 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1966 Mat C; 1967 PetscErrorCode ierr; 1968 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1969 MatScalar *array = a->a; 1970 1971 PetscFunctionBegin; 1972 if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1973 1974 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1975 ierr = PetscCalloc1(1+A->cmap->n,&col);CHKERRQ(ierr); 1976 1977 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1978 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1979 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1980 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 1981 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1982 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1983 ierr = PetscFree(col);CHKERRQ(ierr); 1984 } else { 1985 C = *B; 1986 } 1987 1988 for (i=0; i<m; i++) { 1989 len = ai[i+1]-ai[i]; 1990 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1991 array += len; 1992 aj += len; 1993 } 1994 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1995 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1996 1997 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1998 *B = C; 1999 } else { 2000 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 2001 } 2002 PetscFunctionReturn(0); 2003 } 2004 2005 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2006 { 2007 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2008 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2009 MatScalar *va,*vb; 2010 PetscErrorCode ierr; 2011 PetscInt ma,na,mb,nb, i; 2012 2013 PetscFunctionBegin; 2014 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2015 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2016 if (ma!=nb || na!=mb) { 2017 *f = PETSC_FALSE; 2018 PetscFunctionReturn(0); 2019 } 2020 aii = aij->i; bii = bij->i; 2021 adx = aij->j; bdx = bij->j; 2022 va = aij->a; vb = bij->a; 2023 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2024 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2025 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2026 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2027 2028 *f = PETSC_TRUE; 2029 for (i=0; i<ma; i++) { 2030 while (aptr[i]<aii[i+1]) { 2031 PetscInt idc,idr; 2032 PetscScalar vc,vr; 2033 /* column/row index/value */ 2034 idc = adx[aptr[i]]; 2035 idr = bdx[bptr[idc]]; 2036 vc = va[aptr[i]]; 2037 vr = vb[bptr[idc]]; 2038 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2039 *f = PETSC_FALSE; 2040 goto done; 2041 } else { 2042 aptr[i]++; 2043 if (B || i!=idc) bptr[idc]++; 2044 } 2045 } 2046 } 2047 done: 2048 ierr = PetscFree(aptr);CHKERRQ(ierr); 2049 ierr = PetscFree(bptr);CHKERRQ(ierr); 2050 PetscFunctionReturn(0); 2051 } 2052 2053 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2054 { 2055 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2056 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2057 MatScalar *va,*vb; 2058 PetscErrorCode ierr; 2059 PetscInt ma,na,mb,nb, i; 2060 2061 PetscFunctionBegin; 2062 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2063 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2064 if (ma!=nb || na!=mb) { 2065 *f = PETSC_FALSE; 2066 PetscFunctionReturn(0); 2067 } 2068 aii = aij->i; bii = bij->i; 2069 adx = aij->j; bdx = bij->j; 2070 va = aij->a; vb = bij->a; 2071 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2072 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2073 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2074 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2075 2076 *f = PETSC_TRUE; 2077 for (i=0; i<ma; i++) { 2078 while (aptr[i]<aii[i+1]) { 2079 PetscInt idc,idr; 2080 PetscScalar vc,vr; 2081 /* column/row index/value */ 2082 idc = adx[aptr[i]]; 2083 idr = bdx[bptr[idc]]; 2084 vc = va[aptr[i]]; 2085 vr = vb[bptr[idc]]; 2086 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2087 *f = PETSC_FALSE; 2088 goto done; 2089 } else { 2090 aptr[i]++; 2091 if (B || i!=idc) bptr[idc]++; 2092 } 2093 } 2094 } 2095 done: 2096 ierr = PetscFree(aptr);CHKERRQ(ierr); 2097 ierr = PetscFree(bptr);CHKERRQ(ierr); 2098 PetscFunctionReturn(0); 2099 } 2100 2101 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2102 { 2103 PetscErrorCode ierr; 2104 2105 PetscFunctionBegin; 2106 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2107 PetscFunctionReturn(0); 2108 } 2109 2110 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2111 { 2112 PetscErrorCode ierr; 2113 2114 PetscFunctionBegin; 2115 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2116 PetscFunctionReturn(0); 2117 } 2118 2119 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2120 { 2121 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2122 PetscScalar *l,*r,x; 2123 MatScalar *v; 2124 PetscErrorCode ierr; 2125 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 2126 2127 PetscFunctionBegin; 2128 if (ll) { 2129 /* The local size is used so that VecMPI can be passed to this routine 2130 by MatDiagonalScale_MPIAIJ */ 2131 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2132 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2133 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 2134 v = a->a; 2135 for (i=0; i<m; i++) { 2136 x = l[i]; 2137 M = a->i[i+1] - a->i[i]; 2138 for (j=0; j<M; j++) (*v++) *= x; 2139 } 2140 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 2141 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2142 } 2143 if (rr) { 2144 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2145 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2146 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 2147 v = a->a; jj = a->j; 2148 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2149 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 2150 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2151 } 2152 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2153 PetscFunctionReturn(0); 2154 } 2155 2156 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2157 { 2158 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2159 PetscErrorCode ierr; 2160 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2161 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2162 const PetscInt *irow,*icol; 2163 PetscInt nrows,ncols; 2164 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2165 MatScalar *a_new,*mat_a; 2166 Mat C; 2167 PetscBool stride; 2168 2169 PetscFunctionBegin; 2170 2171 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2172 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2173 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2174 2175 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2176 if (stride) { 2177 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2178 } else { 2179 first = 0; 2180 step = 0; 2181 } 2182 if (stride && step == 1) { 2183 /* special case of contiguous rows */ 2184 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2185 /* loop over new rows determining lens and starting points */ 2186 for (i=0; i<nrows; i++) { 2187 kstart = ai[irow[i]]; 2188 kend = kstart + ailen[irow[i]]; 2189 starts[i] = kstart; 2190 for (k=kstart; k<kend; k++) { 2191 if (aj[k] >= first) { 2192 starts[i] = k; 2193 break; 2194 } 2195 } 2196 sum = 0; 2197 while (k < kend) { 2198 if (aj[k++] >= first+ncols) break; 2199 sum++; 2200 } 2201 lens[i] = sum; 2202 } 2203 /* create submatrix */ 2204 if (scall == MAT_REUSE_MATRIX) { 2205 PetscInt n_cols,n_rows; 2206 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2207 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2208 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2209 C = *B; 2210 } else { 2211 PetscInt rbs,cbs; 2212 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2213 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2214 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2215 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2216 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2217 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2218 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2219 } 2220 c = (Mat_SeqAIJ*)C->data; 2221 2222 /* loop over rows inserting into submatrix */ 2223 a_new = c->a; 2224 j_new = c->j; 2225 i_new = c->i; 2226 2227 for (i=0; i<nrows; i++) { 2228 ii = starts[i]; 2229 lensi = lens[i]; 2230 for (k=0; k<lensi; k++) { 2231 *j_new++ = aj[ii+k] - first; 2232 } 2233 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2234 a_new += lensi; 2235 i_new[i+1] = i_new[i] + lensi; 2236 c->ilen[i] = lensi; 2237 } 2238 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2239 } else { 2240 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2241 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2242 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2243 for (i=0; i<ncols; i++) { 2244 #if defined(PETSC_USE_DEBUG) 2245 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); 2246 #endif 2247 smap[icol[i]] = i+1; 2248 } 2249 2250 /* determine lens of each row */ 2251 for (i=0; i<nrows; i++) { 2252 kstart = ai[irow[i]]; 2253 kend = kstart + a->ilen[irow[i]]; 2254 lens[i] = 0; 2255 for (k=kstart; k<kend; k++) { 2256 if (smap[aj[k]]) { 2257 lens[i]++; 2258 } 2259 } 2260 } 2261 /* Create and fill new matrix */ 2262 if (scall == MAT_REUSE_MATRIX) { 2263 PetscBool equal; 2264 2265 c = (Mat_SeqAIJ*)((*B)->data); 2266 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2267 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2268 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2269 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2270 C = *B; 2271 } else { 2272 PetscInt rbs,cbs; 2273 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2274 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2275 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2276 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2277 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2278 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2279 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2280 } 2281 c = (Mat_SeqAIJ*)(C->data); 2282 for (i=0; i<nrows; i++) { 2283 row = irow[i]; 2284 kstart = ai[row]; 2285 kend = kstart + a->ilen[row]; 2286 mat_i = c->i[i]; 2287 mat_j = c->j + mat_i; 2288 mat_a = c->a + mat_i; 2289 mat_ilen = c->ilen + i; 2290 for (k=kstart; k<kend; k++) { 2291 if ((tcol=smap[a->j[k]])) { 2292 *mat_j++ = tcol - 1; 2293 *mat_a++ = a->a[k]; 2294 (*mat_ilen)++; 2295 2296 } 2297 } 2298 } 2299 /* Free work space */ 2300 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2301 ierr = PetscFree(smap);CHKERRQ(ierr); 2302 ierr = PetscFree(lens);CHKERRQ(ierr); 2303 /* sort */ 2304 for (i = 0; i < nrows; i++) { 2305 PetscInt ilen; 2306 2307 mat_i = c->i[i]; 2308 mat_j = c->j + mat_i; 2309 mat_a = c->a + mat_i; 2310 ilen = c->ilen[i]; 2311 ierr = PetscSortIntWithMatScalarArray(ilen,mat_j,mat_a);CHKERRQ(ierr); 2312 } 2313 } 2314 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2315 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2316 2317 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2318 *B = C; 2319 PetscFunctionReturn(0); 2320 } 2321 2322 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2323 { 2324 PetscErrorCode ierr; 2325 Mat B; 2326 2327 PetscFunctionBegin; 2328 if (scall == MAT_INITIAL_MATRIX) { 2329 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2330 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2331 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2332 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2333 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2334 *subMat = B; 2335 } else { 2336 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2337 } 2338 PetscFunctionReturn(0); 2339 } 2340 2341 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2342 { 2343 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2344 PetscErrorCode ierr; 2345 Mat outA; 2346 PetscBool row_identity,col_identity; 2347 2348 PetscFunctionBegin; 2349 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2350 2351 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2352 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2353 2354 outA = inA; 2355 outA->factortype = MAT_FACTOR_LU; 2356 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 2357 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 2358 2359 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2360 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2361 2362 a->row = row; 2363 2364 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2365 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2366 2367 a->col = col; 2368 2369 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2370 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2371 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2372 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2373 2374 if (!a->solve_work) { /* this matrix may have been factored before */ 2375 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2376 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2377 } 2378 2379 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2380 if (row_identity && col_identity) { 2381 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2382 } else { 2383 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2384 } 2385 PetscFunctionReturn(0); 2386 } 2387 2388 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2389 { 2390 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2391 PetscScalar oalpha = alpha; 2392 PetscErrorCode ierr; 2393 PetscBLASInt one = 1,bnz; 2394 2395 PetscFunctionBegin; 2396 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2397 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2398 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2399 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2400 PetscFunctionReturn(0); 2401 } 2402 2403 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2404 { 2405 PetscErrorCode ierr; 2406 PetscInt i; 2407 2408 PetscFunctionBegin; 2409 if (scall == MAT_INITIAL_MATRIX) { 2410 ierr = PetscMalloc1(n+1,B);CHKERRQ(ierr); 2411 } 2412 2413 for (i=0; i<n; i++) { 2414 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2415 } 2416 PetscFunctionReturn(0); 2417 } 2418 2419 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2420 { 2421 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2422 PetscErrorCode ierr; 2423 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2424 const PetscInt *idx; 2425 PetscInt start,end,*ai,*aj; 2426 PetscBT table; 2427 2428 PetscFunctionBegin; 2429 m = A->rmap->n; 2430 ai = a->i; 2431 aj = a->j; 2432 2433 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2434 2435 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2436 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2437 2438 for (i=0; i<is_max; i++) { 2439 /* Initialize the two local arrays */ 2440 isz = 0; 2441 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2442 2443 /* Extract the indices, assume there can be duplicate entries */ 2444 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2445 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2446 2447 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2448 for (j=0; j<n; ++j) { 2449 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2450 } 2451 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2452 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2453 2454 k = 0; 2455 for (j=0; j<ov; j++) { /* for each overlap */ 2456 n = isz; 2457 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2458 row = nidx[k]; 2459 start = ai[row]; 2460 end = ai[row+1]; 2461 for (l = start; l<end; l++) { 2462 val = aj[l]; 2463 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2464 } 2465 } 2466 } 2467 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2468 } 2469 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2470 ierr = PetscFree(nidx);CHKERRQ(ierr); 2471 PetscFunctionReturn(0); 2472 } 2473 2474 /* -------------------------------------------------------------- */ 2475 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2476 { 2477 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2478 PetscErrorCode ierr; 2479 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2480 const PetscInt *row,*col; 2481 PetscInt *cnew,j,*lens; 2482 IS icolp,irowp; 2483 PetscInt *cwork = NULL; 2484 PetscScalar *vwork = NULL; 2485 2486 PetscFunctionBegin; 2487 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2488 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2489 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2490 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2491 2492 /* determine lengths of permuted rows */ 2493 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2494 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2495 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2496 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2497 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2498 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2499 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2500 ierr = PetscFree(lens);CHKERRQ(ierr); 2501 2502 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2503 for (i=0; i<m; i++) { 2504 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2505 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2506 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2507 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2508 } 2509 ierr = PetscFree(cnew);CHKERRQ(ierr); 2510 2511 (*B)->assembled = PETSC_FALSE; 2512 2513 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2514 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2515 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2516 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2517 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2518 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2519 PetscFunctionReturn(0); 2520 } 2521 2522 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2523 { 2524 PetscErrorCode ierr; 2525 2526 PetscFunctionBegin; 2527 /* If the two matrices have the same copy implementation, use fast copy. */ 2528 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2529 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2530 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2531 2532 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"); 2533 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2534 } else { 2535 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2536 } 2537 PetscFunctionReturn(0); 2538 } 2539 2540 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2541 { 2542 PetscErrorCode ierr; 2543 2544 PetscFunctionBegin; 2545 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2546 PetscFunctionReturn(0); 2547 } 2548 2549 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2550 { 2551 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2552 2553 PetscFunctionBegin; 2554 *array = a->a; 2555 PetscFunctionReturn(0); 2556 } 2557 2558 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2559 { 2560 PetscFunctionBegin; 2561 PetscFunctionReturn(0); 2562 } 2563 2564 /* 2565 Computes the number of nonzeros per row needed for preallocation when X and Y 2566 have different nonzero structure. 2567 */ 2568 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2569 { 2570 PetscInt i,j,k,nzx,nzy; 2571 2572 PetscFunctionBegin; 2573 /* Set the number of nonzeros in the new matrix */ 2574 for (i=0; i<m; i++) { 2575 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2576 nzx = xi[i+1] - xi[i]; 2577 nzy = yi[i+1] - yi[i]; 2578 nnz[i] = 0; 2579 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2580 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2581 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2582 nnz[i]++; 2583 } 2584 for (; k<nzy; k++) nnz[i]++; 2585 } 2586 PetscFunctionReturn(0); 2587 } 2588 2589 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2590 { 2591 PetscInt m = Y->rmap->N; 2592 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2593 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2594 PetscErrorCode ierr; 2595 2596 PetscFunctionBegin; 2597 /* Set the number of nonzeros in the new matrix */ 2598 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2599 PetscFunctionReturn(0); 2600 } 2601 2602 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2603 { 2604 PetscErrorCode ierr; 2605 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2606 PetscBLASInt one=1,bnz; 2607 2608 PetscFunctionBegin; 2609 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2610 if (str == SAME_NONZERO_PATTERN) { 2611 PetscScalar alpha = a; 2612 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2613 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2614 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2615 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2616 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2617 } else { 2618 Mat B; 2619 PetscInt *nnz; 2620 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2621 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2622 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2623 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2624 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2625 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2626 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2627 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2628 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2629 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2630 ierr = PetscFree(nnz);CHKERRQ(ierr); 2631 } 2632 PetscFunctionReturn(0); 2633 } 2634 2635 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2636 { 2637 #if defined(PETSC_USE_COMPLEX) 2638 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2639 PetscInt i,nz; 2640 PetscScalar *a; 2641 2642 PetscFunctionBegin; 2643 nz = aij->nz; 2644 a = aij->a; 2645 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2646 #else 2647 PetscFunctionBegin; 2648 #endif 2649 PetscFunctionReturn(0); 2650 } 2651 2652 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2653 { 2654 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2655 PetscErrorCode ierr; 2656 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2657 PetscReal atmp; 2658 PetscScalar *x; 2659 MatScalar *aa; 2660 2661 PetscFunctionBegin; 2662 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2663 aa = a->a; 2664 ai = a->i; 2665 aj = a->j; 2666 2667 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2668 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2669 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2670 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2671 for (i=0; i<m; i++) { 2672 ncols = ai[1] - ai[0]; ai++; 2673 x[i] = 0.0; 2674 for (j=0; j<ncols; j++) { 2675 atmp = PetscAbsScalar(*aa); 2676 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2677 aa++; aj++; 2678 } 2679 } 2680 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2681 PetscFunctionReturn(0); 2682 } 2683 2684 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2685 { 2686 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2687 PetscErrorCode ierr; 2688 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2689 PetscScalar *x; 2690 MatScalar *aa; 2691 2692 PetscFunctionBegin; 2693 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2694 aa = a->a; 2695 ai = a->i; 2696 aj = a->j; 2697 2698 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2699 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2700 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2701 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2702 for (i=0; i<m; i++) { 2703 ncols = ai[1] - ai[0]; ai++; 2704 if (ncols == A->cmap->n) { /* row is dense */ 2705 x[i] = *aa; if (idx) idx[i] = 0; 2706 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2707 x[i] = 0.0; 2708 if (idx) { 2709 idx[i] = 0; /* in case ncols is zero */ 2710 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2711 if (aj[j] > j) { 2712 idx[i] = j; 2713 break; 2714 } 2715 } 2716 } 2717 } 2718 for (j=0; j<ncols; j++) { 2719 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2720 aa++; aj++; 2721 } 2722 } 2723 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2724 PetscFunctionReturn(0); 2725 } 2726 2727 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2728 { 2729 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2730 PetscErrorCode ierr; 2731 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2732 PetscReal atmp; 2733 PetscScalar *x; 2734 MatScalar *aa; 2735 2736 PetscFunctionBegin; 2737 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2738 aa = a->a; 2739 ai = a->i; 2740 aj = a->j; 2741 2742 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2743 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2744 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2745 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); 2746 for (i=0; i<m; i++) { 2747 ncols = ai[1] - ai[0]; ai++; 2748 if (ncols) { 2749 /* Get first nonzero */ 2750 for (j = 0; j < ncols; j++) { 2751 atmp = PetscAbsScalar(aa[j]); 2752 if (atmp > 1.0e-12) { 2753 x[i] = atmp; 2754 if (idx) idx[i] = aj[j]; 2755 break; 2756 } 2757 } 2758 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2759 } else { 2760 x[i] = 0.0; if (idx) idx[i] = 0; 2761 } 2762 for (j = 0; j < ncols; j++) { 2763 atmp = PetscAbsScalar(*aa); 2764 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2765 aa++; aj++; 2766 } 2767 } 2768 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2769 PetscFunctionReturn(0); 2770 } 2771 2772 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2773 { 2774 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2775 PetscErrorCode ierr; 2776 PetscInt i,j,m = A->rmap->n,ncols,n; 2777 const PetscInt *ai,*aj; 2778 PetscScalar *x; 2779 const MatScalar *aa; 2780 2781 PetscFunctionBegin; 2782 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2783 aa = a->a; 2784 ai = a->i; 2785 aj = a->j; 2786 2787 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2788 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2789 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2790 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2791 for (i=0; i<m; i++) { 2792 ncols = ai[1] - ai[0]; ai++; 2793 if (ncols == A->cmap->n) { /* row is dense */ 2794 x[i] = *aa; if (idx) idx[i] = 0; 2795 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2796 x[i] = 0.0; 2797 if (idx) { /* find first implicit 0.0 in the row */ 2798 idx[i] = 0; /* in case ncols is zero */ 2799 for (j=0; j<ncols; j++) { 2800 if (aj[j] > j) { 2801 idx[i] = j; 2802 break; 2803 } 2804 } 2805 } 2806 } 2807 for (j=0; j<ncols; j++) { 2808 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2809 aa++; aj++; 2810 } 2811 } 2812 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2813 PetscFunctionReturn(0); 2814 } 2815 2816 #include <petscblaslapack.h> 2817 #include <petsc/private/kernels/blockinvert.h> 2818 2819 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 2820 { 2821 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 2822 PetscErrorCode ierr; 2823 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 2824 MatScalar *diag,work[25],*v_work; 2825 PetscReal shift = 0.0; 2826 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 2827 2828 PetscFunctionBegin; 2829 allowzeropivot = PetscNot(A->erroriffailure); 2830 if (a->ibdiagvalid) { 2831 if (values) *values = a->ibdiag; 2832 PetscFunctionReturn(0); 2833 } 2834 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 2835 if (!a->ibdiag) { 2836 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 2837 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2838 } 2839 diag = a->ibdiag; 2840 if (values) *values = a->ibdiag; 2841 /* factor and invert each block */ 2842 switch (bs) { 2843 case 1: 2844 for (i=0; i<mbs; i++) { 2845 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 2846 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 2847 if (allowzeropivot) { 2848 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2849 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 2850 A->factorerror_zeropivot_row = i; 2851 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 2852 } 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); 2853 } 2854 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 2855 } 2856 break; 2857 case 2: 2858 for (i=0; i<mbs; i++) { 2859 ij[0] = 2*i; ij[1] = 2*i + 1; 2860 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 2861 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2862 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2863 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 2864 diag += 4; 2865 } 2866 break; 2867 case 3: 2868 for (i=0; i<mbs; i++) { 2869 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 2870 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 2871 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2872 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2873 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 2874 diag += 9; 2875 } 2876 break; 2877 case 4: 2878 for (i=0; i<mbs; i++) { 2879 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 2880 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 2881 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2882 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2883 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 2884 diag += 16; 2885 } 2886 break; 2887 case 5: 2888 for (i=0; i<mbs; i++) { 2889 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 2890 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 2891 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2892 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2893 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 2894 diag += 25; 2895 } 2896 break; 2897 case 6: 2898 for (i=0; i<mbs; i++) { 2899 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; 2900 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 2901 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2902 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2903 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 2904 diag += 36; 2905 } 2906 break; 2907 case 7: 2908 for (i=0; i<mbs; i++) { 2909 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; 2910 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 2911 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2912 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2913 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 2914 diag += 49; 2915 } 2916 break; 2917 default: 2918 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 2919 for (i=0; i<mbs; i++) { 2920 for (j=0; j<bs; j++) { 2921 IJ[j] = bs*i + j; 2922 } 2923 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 2924 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2925 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2926 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 2927 diag += bs2; 2928 } 2929 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 2930 } 2931 a->ibdiagvalid = PETSC_TRUE; 2932 PetscFunctionReturn(0); 2933 } 2934 2935 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 2936 { 2937 PetscErrorCode ierr; 2938 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 2939 PetscScalar a; 2940 PetscInt m,n,i,j,col; 2941 2942 PetscFunctionBegin; 2943 if (!x->assembled) { 2944 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2945 for (i=0; i<m; i++) { 2946 for (j=0; j<aij->imax[i]; j++) { 2947 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 2948 col = (PetscInt)(n*PetscRealPart(a)); 2949 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 2950 } 2951 } 2952 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 2953 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2954 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2955 PetscFunctionReturn(0); 2956 } 2957 2958 PetscErrorCode MatShift_SeqAIJ(Mat Y,PetscScalar a) 2959 { 2960 PetscErrorCode ierr; 2961 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)Y->data; 2962 2963 PetscFunctionBegin; 2964 if (!Y->preallocated || !aij->nz) { 2965 ierr = MatSeqAIJSetPreallocation(Y,1,NULL);CHKERRQ(ierr); 2966 } 2967 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 2968 PetscFunctionReturn(0); 2969 } 2970 2971 /* -------------------------------------------------------------------*/ 2972 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 2973 MatGetRow_SeqAIJ, 2974 MatRestoreRow_SeqAIJ, 2975 MatMult_SeqAIJ, 2976 /* 4*/ MatMultAdd_SeqAIJ, 2977 MatMultTranspose_SeqAIJ, 2978 MatMultTransposeAdd_SeqAIJ, 2979 0, 2980 0, 2981 0, 2982 /* 10*/ 0, 2983 MatLUFactor_SeqAIJ, 2984 0, 2985 MatSOR_SeqAIJ, 2986 MatTranspose_SeqAIJ, 2987 /*1 5*/ MatGetInfo_SeqAIJ, 2988 MatEqual_SeqAIJ, 2989 MatGetDiagonal_SeqAIJ, 2990 MatDiagonalScale_SeqAIJ, 2991 MatNorm_SeqAIJ, 2992 /* 20*/ 0, 2993 MatAssemblyEnd_SeqAIJ, 2994 MatSetOption_SeqAIJ, 2995 MatZeroEntries_SeqAIJ, 2996 /* 24*/ MatZeroRows_SeqAIJ, 2997 0, 2998 0, 2999 0, 3000 0, 3001 /* 29*/ MatSetUp_SeqAIJ, 3002 0, 3003 0, 3004 0, 3005 0, 3006 /* 34*/ MatDuplicate_SeqAIJ, 3007 0, 3008 0, 3009 MatILUFactor_SeqAIJ, 3010 0, 3011 /* 39*/ MatAXPY_SeqAIJ, 3012 MatGetSubMatrices_SeqAIJ, 3013 MatIncreaseOverlap_SeqAIJ, 3014 MatGetValues_SeqAIJ, 3015 MatCopy_SeqAIJ, 3016 /* 44*/ MatGetRowMax_SeqAIJ, 3017 MatScale_SeqAIJ, 3018 MatShift_SeqAIJ, 3019 MatDiagonalSet_SeqAIJ, 3020 MatZeroRowsColumns_SeqAIJ, 3021 /* 49*/ MatSetRandom_SeqAIJ, 3022 MatGetRowIJ_SeqAIJ, 3023 MatRestoreRowIJ_SeqAIJ, 3024 MatGetColumnIJ_SeqAIJ, 3025 MatRestoreColumnIJ_SeqAIJ, 3026 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3027 0, 3028 0, 3029 MatPermute_SeqAIJ, 3030 0, 3031 /* 59*/ 0, 3032 MatDestroy_SeqAIJ, 3033 MatView_SeqAIJ, 3034 0, 3035 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3036 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3037 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3038 0, 3039 0, 3040 0, 3041 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3042 MatGetRowMinAbs_SeqAIJ, 3043 0, 3044 0, 3045 0, 3046 /* 74*/ 0, 3047 MatFDColoringApply_AIJ, 3048 0, 3049 0, 3050 0, 3051 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3052 0, 3053 0, 3054 0, 3055 MatLoad_SeqAIJ, 3056 /* 84*/ MatIsSymmetric_SeqAIJ, 3057 MatIsHermitian_SeqAIJ, 3058 0, 3059 0, 3060 0, 3061 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3062 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3063 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3064 MatPtAP_SeqAIJ_SeqAIJ, 3065 MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy, 3066 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3067 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3068 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3069 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3070 0, 3071 /* 99*/ 0, 3072 0, 3073 0, 3074 MatConjugate_SeqAIJ, 3075 0, 3076 /*104*/ MatSetValuesRow_SeqAIJ, 3077 MatRealPart_SeqAIJ, 3078 MatImaginaryPart_SeqAIJ, 3079 0, 3080 0, 3081 /*109*/ MatMatSolve_SeqAIJ, 3082 0, 3083 MatGetRowMin_SeqAIJ, 3084 0, 3085 MatMissingDiagonal_SeqAIJ, 3086 /*114*/ 0, 3087 0, 3088 0, 3089 0, 3090 0, 3091 /*119*/ 0, 3092 0, 3093 0, 3094 0, 3095 MatGetMultiProcBlock_SeqAIJ, 3096 /*124*/ MatFindNonzeroRows_SeqAIJ, 3097 MatGetColumnNorms_SeqAIJ, 3098 MatInvertBlockDiagonal_SeqAIJ, 3099 0, 3100 0, 3101 /*129*/ 0, 3102 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3103 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3104 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3105 MatTransposeColoringCreate_SeqAIJ, 3106 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3107 MatTransColoringApplyDenToSp_SeqAIJ, 3108 MatRARt_SeqAIJ_SeqAIJ, 3109 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3110 MatRARtNumeric_SeqAIJ_SeqAIJ, 3111 /*139*/0, 3112 0, 3113 0, 3114 MatFDColoringSetUp_SeqXAIJ, 3115 MatFindOffBlockDiagonalEntries_SeqAIJ, 3116 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ 3117 }; 3118 3119 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3120 { 3121 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3122 PetscInt i,nz,n; 3123 3124 PetscFunctionBegin; 3125 nz = aij->maxnz; 3126 n = mat->rmap->n; 3127 for (i=0; i<nz; i++) { 3128 aij->j[i] = indices[i]; 3129 } 3130 aij->nz = nz; 3131 for (i=0; i<n; i++) { 3132 aij->ilen[i] = aij->imax[i]; 3133 } 3134 PetscFunctionReturn(0); 3135 } 3136 3137 /*@ 3138 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3139 in the matrix. 3140 3141 Input Parameters: 3142 + mat - the SeqAIJ matrix 3143 - indices - the column indices 3144 3145 Level: advanced 3146 3147 Notes: 3148 This can be called if you have precomputed the nonzero structure of the 3149 matrix and want to provide it to the matrix object to improve the performance 3150 of the MatSetValues() operation. 3151 3152 You MUST have set the correct numbers of nonzeros per row in the call to 3153 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3154 3155 MUST be called before any calls to MatSetValues(); 3156 3157 The indices should start with zero, not one. 3158 3159 @*/ 3160 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3161 { 3162 PetscErrorCode ierr; 3163 3164 PetscFunctionBegin; 3165 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3166 PetscValidPointer(indices,2); 3167 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3168 PetscFunctionReturn(0); 3169 } 3170 3171 /* ----------------------------------------------------------------------------------------*/ 3172 3173 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3174 { 3175 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3176 PetscErrorCode ierr; 3177 size_t nz = aij->i[mat->rmap->n]; 3178 3179 PetscFunctionBegin; 3180 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3181 3182 /* allocate space for values if not already there */ 3183 if (!aij->saved_values) { 3184 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3185 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3186 } 3187 3188 /* copy values over */ 3189 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3190 PetscFunctionReturn(0); 3191 } 3192 3193 /*@ 3194 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3195 example, reuse of the linear part of a Jacobian, while recomputing the 3196 nonlinear portion. 3197 3198 Collect on Mat 3199 3200 Input Parameters: 3201 . mat - the matrix (currently only AIJ matrices support this option) 3202 3203 Level: advanced 3204 3205 Common Usage, with SNESSolve(): 3206 $ Create Jacobian matrix 3207 $ Set linear terms into matrix 3208 $ Apply boundary conditions to matrix, at this time matrix must have 3209 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3210 $ boundary conditions again will not change the nonzero structure 3211 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3212 $ ierr = MatStoreValues(mat); 3213 $ Call SNESSetJacobian() with matrix 3214 $ In your Jacobian routine 3215 $ ierr = MatRetrieveValues(mat); 3216 $ Set nonlinear terms in matrix 3217 3218 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3219 $ // build linear portion of Jacobian 3220 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3221 $ ierr = MatStoreValues(mat); 3222 $ loop over nonlinear iterations 3223 $ ierr = MatRetrieveValues(mat); 3224 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3225 $ // call MatAssemblyBegin/End() on matrix 3226 $ Solve linear system with Jacobian 3227 $ endloop 3228 3229 Notes: 3230 Matrix must already be assemblied before calling this routine 3231 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3232 calling this routine. 3233 3234 When this is called multiple times it overwrites the previous set of stored values 3235 and does not allocated additional space. 3236 3237 .seealso: MatRetrieveValues() 3238 3239 @*/ 3240 PetscErrorCode MatStoreValues(Mat mat) 3241 { 3242 PetscErrorCode ierr; 3243 3244 PetscFunctionBegin; 3245 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3246 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3247 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3248 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3249 PetscFunctionReturn(0); 3250 } 3251 3252 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3253 { 3254 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3255 PetscErrorCode ierr; 3256 PetscInt nz = aij->i[mat->rmap->n]; 3257 3258 PetscFunctionBegin; 3259 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3260 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3261 /* copy values over */ 3262 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3263 PetscFunctionReturn(0); 3264 } 3265 3266 /*@ 3267 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3268 example, reuse of the linear part of a Jacobian, while recomputing the 3269 nonlinear portion. 3270 3271 Collect on Mat 3272 3273 Input Parameters: 3274 . mat - the matrix (currently only AIJ matrices support this option) 3275 3276 Level: advanced 3277 3278 .seealso: MatStoreValues() 3279 3280 @*/ 3281 PetscErrorCode MatRetrieveValues(Mat mat) 3282 { 3283 PetscErrorCode ierr; 3284 3285 PetscFunctionBegin; 3286 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3287 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3288 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3289 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3290 PetscFunctionReturn(0); 3291 } 3292 3293 3294 /* --------------------------------------------------------------------------------*/ 3295 /*@C 3296 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3297 (the default parallel PETSc format). For good matrix assembly performance 3298 the user should preallocate the matrix storage by setting the parameter nz 3299 (or the array nnz). By setting these parameters accurately, performance 3300 during matrix assembly can be increased by more than a factor of 50. 3301 3302 Collective on MPI_Comm 3303 3304 Input Parameters: 3305 + comm - MPI communicator, set to PETSC_COMM_SELF 3306 . m - number of rows 3307 . n - number of columns 3308 . nz - number of nonzeros per row (same for all rows) 3309 - nnz - array containing the number of nonzeros in the various rows 3310 (possibly different for each row) or NULL 3311 3312 Output Parameter: 3313 . A - the matrix 3314 3315 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3316 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3317 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3318 3319 Notes: 3320 If nnz is given then nz is ignored 3321 3322 The AIJ format (also called the Yale sparse matrix format or 3323 compressed row storage), is fully compatible with standard Fortran 77 3324 storage. That is, the stored row and column indices can begin at 3325 either one (as in Fortran) or zero. See the users' manual for details. 3326 3327 Specify the preallocated storage with either nz or nnz (not both). 3328 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3329 allocation. For large problems you MUST preallocate memory or you 3330 will get TERRIBLE performance, see the users' manual chapter on matrices. 3331 3332 By default, this format uses inodes (identical nodes) when possible, to 3333 improve numerical efficiency of matrix-vector products and solves. We 3334 search for consecutive rows with the same nonzero structure, thereby 3335 reusing matrix information to achieve increased efficiency. 3336 3337 Options Database Keys: 3338 + -mat_no_inode - Do not use inodes 3339 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3340 3341 Level: intermediate 3342 3343 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3344 3345 @*/ 3346 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3347 { 3348 PetscErrorCode ierr; 3349 3350 PetscFunctionBegin; 3351 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3352 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3353 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3354 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3355 PetscFunctionReturn(0); 3356 } 3357 3358 /*@C 3359 MatSeqAIJSetPreallocation - For good matrix assembly performance 3360 the user should preallocate the matrix storage by setting the parameter nz 3361 (or the array nnz). By setting these parameters accurately, performance 3362 during matrix assembly can be increased by more than a factor of 50. 3363 3364 Collective on MPI_Comm 3365 3366 Input Parameters: 3367 + B - The matrix 3368 . nz - number of nonzeros per row (same for all rows) 3369 - nnz - array containing the number of nonzeros in the various rows 3370 (possibly different for each row) or NULL 3371 3372 Notes: 3373 If nnz is given then nz is ignored 3374 3375 The AIJ format (also called the Yale sparse matrix format or 3376 compressed row storage), is fully compatible with standard Fortran 77 3377 storage. That is, the stored row and column indices can begin at 3378 either one (as in Fortran) or zero. See the users' manual for details. 3379 3380 Specify the preallocated storage with either nz or nnz (not both). 3381 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3382 allocation. For large problems you MUST preallocate memory or you 3383 will get TERRIBLE performance, see the users' manual chapter on matrices. 3384 3385 You can call MatGetInfo() to get information on how effective the preallocation was; 3386 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3387 You can also run with the option -info and look for messages with the string 3388 malloc in them to see if additional memory allocation was needed. 3389 3390 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3391 entries or columns indices 3392 3393 By default, this format uses inodes (identical nodes) when possible, to 3394 improve numerical efficiency of matrix-vector products and solves. We 3395 search for consecutive rows with the same nonzero structure, thereby 3396 reusing matrix information to achieve increased efficiency. 3397 3398 Options Database Keys: 3399 + -mat_no_inode - Do not use inodes 3400 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3401 - -mat_aij_oneindex - Internally use indexing starting at 1 3402 rather than 0. Note that when calling MatSetValues(), 3403 the user still MUST index entries starting at 0! 3404 3405 Level: intermediate 3406 3407 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3408 3409 @*/ 3410 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3411 { 3412 PetscErrorCode ierr; 3413 3414 PetscFunctionBegin; 3415 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3416 PetscValidType(B,1); 3417 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3418 PetscFunctionReturn(0); 3419 } 3420 3421 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3422 { 3423 Mat_SeqAIJ *b; 3424 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3425 PetscErrorCode ierr; 3426 PetscInt i; 3427 3428 PetscFunctionBegin; 3429 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3430 if (nz == MAT_SKIP_ALLOCATION) { 3431 skipallocation = PETSC_TRUE; 3432 nz = 0; 3433 } 3434 3435 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3436 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3437 3438 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3439 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3440 if (nnz) { 3441 for (i=0; i<B->rmap->n; i++) { 3442 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]); 3443 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); 3444 } 3445 } 3446 3447 B->preallocated = PETSC_TRUE; 3448 3449 b = (Mat_SeqAIJ*)B->data; 3450 3451 if (!skipallocation) { 3452 if (!b->imax) { 3453 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3454 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3455 } 3456 if (!nnz) { 3457 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3458 else if (nz < 0) nz = 1; 3459 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3460 nz = nz*B->rmap->n; 3461 } else { 3462 nz = 0; 3463 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3464 } 3465 /* b->ilen will count nonzeros in each row so far. */ 3466 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3467 3468 /* allocate the matrix space */ 3469 /* FIXME: should B's old memory be unlogged? */ 3470 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3471 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3472 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3473 b->i[0] = 0; 3474 for (i=1; i<B->rmap->n+1; i++) { 3475 b->i[i] = b->i[i-1] + b->imax[i-1]; 3476 } 3477 b->singlemalloc = PETSC_TRUE; 3478 b->free_a = PETSC_TRUE; 3479 b->free_ij = PETSC_TRUE; 3480 } else { 3481 b->free_a = PETSC_FALSE; 3482 b->free_ij = PETSC_FALSE; 3483 } 3484 3485 b->nz = 0; 3486 b->maxnz = nz; 3487 B->info.nz_unneeded = (double)b->maxnz; 3488 if (realalloc) { 3489 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3490 } 3491 B->was_assembled = PETSC_FALSE; 3492 B->assembled = PETSC_FALSE; 3493 PetscFunctionReturn(0); 3494 } 3495 3496 /*@ 3497 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3498 3499 Input Parameters: 3500 + B - the matrix 3501 . i - the indices into j for the start of each row (starts with zero) 3502 . j - the column indices for each row (starts with zero) these must be sorted for each row 3503 - v - optional values in the matrix 3504 3505 Level: developer 3506 3507 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3508 3509 .keywords: matrix, aij, compressed row, sparse, sequential 3510 3511 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3512 @*/ 3513 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3514 { 3515 PetscErrorCode ierr; 3516 3517 PetscFunctionBegin; 3518 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3519 PetscValidType(B,1); 3520 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3521 PetscFunctionReturn(0); 3522 } 3523 3524 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3525 { 3526 PetscInt i; 3527 PetscInt m,n; 3528 PetscInt nz; 3529 PetscInt *nnz, nz_max = 0; 3530 PetscScalar *values; 3531 PetscErrorCode ierr; 3532 3533 PetscFunctionBegin; 3534 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3535 3536 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3537 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3538 3539 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3540 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3541 for (i = 0; i < m; i++) { 3542 nz = Ii[i+1]- Ii[i]; 3543 nz_max = PetscMax(nz_max, nz); 3544 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3545 nnz[i] = nz; 3546 } 3547 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3548 ierr = PetscFree(nnz);CHKERRQ(ierr); 3549 3550 if (v) { 3551 values = (PetscScalar*) v; 3552 } else { 3553 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3554 } 3555 3556 for (i = 0; i < m; i++) { 3557 nz = Ii[i+1] - Ii[i]; 3558 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3559 } 3560 3561 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3562 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3563 3564 if (!v) { 3565 ierr = PetscFree(values);CHKERRQ(ierr); 3566 } 3567 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3568 PetscFunctionReturn(0); 3569 } 3570 3571 #include <../src/mat/impls/dense/seq/dense.h> 3572 #include <petsc/private/kernels/petscaxpy.h> 3573 3574 /* 3575 Computes (B'*A')' since computing B*A directly is untenable 3576 3577 n p p 3578 ( ) ( ) ( ) 3579 m ( A ) * n ( B ) = m ( C ) 3580 ( ) ( ) ( ) 3581 3582 */ 3583 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3584 { 3585 PetscErrorCode ierr; 3586 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3587 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3588 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3589 PetscInt i,n,m,q,p; 3590 const PetscInt *ii,*idx; 3591 const PetscScalar *b,*a,*a_q; 3592 PetscScalar *c,*c_q; 3593 3594 PetscFunctionBegin; 3595 m = A->rmap->n; 3596 n = A->cmap->n; 3597 p = B->cmap->n; 3598 a = sub_a->v; 3599 b = sub_b->a; 3600 c = sub_c->v; 3601 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3602 3603 ii = sub_b->i; 3604 idx = sub_b->j; 3605 for (i=0; i<n; i++) { 3606 q = ii[i+1] - ii[i]; 3607 while (q-->0) { 3608 c_q = c + m*(*idx); 3609 a_q = a + m*i; 3610 PetscKernelAXPY(c_q,*b,a_q,m); 3611 idx++; 3612 b++; 3613 } 3614 } 3615 PetscFunctionReturn(0); 3616 } 3617 3618 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3619 { 3620 PetscErrorCode ierr; 3621 PetscInt m=A->rmap->n,n=B->cmap->n; 3622 Mat Cmat; 3623 3624 PetscFunctionBegin; 3625 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); 3626 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3627 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3628 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3629 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3630 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3631 3632 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3633 3634 *C = Cmat; 3635 PetscFunctionReturn(0); 3636 } 3637 3638 /* ----------------------------------------------------------------*/ 3639 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3640 { 3641 PetscErrorCode ierr; 3642 3643 PetscFunctionBegin; 3644 if (scall == MAT_INITIAL_MATRIX) { 3645 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3646 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3647 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3648 } 3649 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3650 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3651 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3652 PetscFunctionReturn(0); 3653 } 3654 3655 3656 /*MC 3657 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3658 based on compressed sparse row format. 3659 3660 Options Database Keys: 3661 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3662 3663 Level: beginner 3664 3665 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3666 M*/ 3667 3668 /*MC 3669 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3670 3671 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3672 and MATMPIAIJ otherwise. As a result, for single process communicators, 3673 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3674 for communicators controlling multiple processes. It is recommended that you call both of 3675 the above preallocation routines for simplicity. 3676 3677 Options Database Keys: 3678 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3679 3680 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3681 enough exist. 3682 3683 Level: beginner 3684 3685 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3686 M*/ 3687 3688 /*MC 3689 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3690 3691 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3692 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3693 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3694 for communicators controlling multiple processes. It is recommended that you call both of 3695 the above preallocation routines for simplicity. 3696 3697 Options Database Keys: 3698 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3699 3700 Level: beginner 3701 3702 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3703 M*/ 3704 3705 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3706 #if defined(PETSC_HAVE_ELEMENTAL) 3707 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3708 #endif 3709 #if defined(PETSC_HAVE_HYPRE) 3710 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 3711 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 3712 #endif 3713 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 3714 3715 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3716 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3717 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3718 #endif 3719 3720 3721 /*@C 3722 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 3723 3724 Not Collective 3725 3726 Input Parameter: 3727 . mat - a MATSEQAIJ matrix 3728 3729 Output Parameter: 3730 . array - pointer to the data 3731 3732 Level: intermediate 3733 3734 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3735 @*/ 3736 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3737 { 3738 PetscErrorCode ierr; 3739 3740 PetscFunctionBegin; 3741 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3742 PetscFunctionReturn(0); 3743 } 3744 3745 /*@C 3746 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 3747 3748 Not Collective 3749 3750 Input Parameter: 3751 . mat - a MATSEQAIJ matrix 3752 3753 Output Parameter: 3754 . nz - the maximum number of nonzeros in any row 3755 3756 Level: intermediate 3757 3758 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3759 @*/ 3760 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 3761 { 3762 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 3763 3764 PetscFunctionBegin; 3765 *nz = aij->rmax; 3766 PetscFunctionReturn(0); 3767 } 3768 3769 /*@C 3770 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 3771 3772 Not Collective 3773 3774 Input Parameters: 3775 . mat - a MATSEQAIJ matrix 3776 . array - pointer to the data 3777 3778 Level: intermediate 3779 3780 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 3781 @*/ 3782 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3783 { 3784 PetscErrorCode ierr; 3785 3786 PetscFunctionBegin; 3787 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3788 PetscFunctionReturn(0); 3789 } 3790 3791 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 3792 { 3793 Mat_SeqAIJ *b; 3794 PetscErrorCode ierr; 3795 PetscMPIInt size; 3796 3797 PetscFunctionBegin; 3798 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3799 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3800 3801 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3802 3803 B->data = (void*)b; 3804 3805 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3806 3807 b->row = 0; 3808 b->col = 0; 3809 b->icol = 0; 3810 b->reallocs = 0; 3811 b->ignorezeroentries = PETSC_FALSE; 3812 b->roworiented = PETSC_TRUE; 3813 b->nonew = 0; 3814 b->diag = 0; 3815 b->solve_work = 0; 3816 B->spptr = 0; 3817 b->saved_values = 0; 3818 b->idiag = 0; 3819 b->mdiag = 0; 3820 b->ssor_work = 0; 3821 b->omega = 1.0; 3822 b->fshift = 0.0; 3823 b->idiagvalid = PETSC_FALSE; 3824 b->ibdiagvalid = PETSC_FALSE; 3825 b->keepnonzeropattern = PETSC_FALSE; 3826 3827 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3828 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 3829 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 3830 3831 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3832 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3833 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3834 #endif 3835 3836 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3837 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3838 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3839 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3840 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3841 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3842 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3843 #if defined(PETSC_HAVE_ELEMENTAL) 3844 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 3845 #endif 3846 #if defined(PETSC_HAVE_HYPRE) 3847 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 3848 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 3849 #endif 3850 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 3851 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3852 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3853 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3854 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3855 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3856 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3857 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3858 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3859 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3860 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3861 PetscFunctionReturn(0); 3862 } 3863 3864 /* 3865 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3866 */ 3867 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3868 { 3869 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3870 PetscErrorCode ierr; 3871 PetscInt i,m = A->rmap->n; 3872 3873 PetscFunctionBegin; 3874 c = (Mat_SeqAIJ*)C->data; 3875 3876 C->factortype = A->factortype; 3877 c->row = 0; 3878 c->col = 0; 3879 c->icol = 0; 3880 c->reallocs = 0; 3881 3882 C->assembled = PETSC_TRUE; 3883 3884 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 3885 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 3886 3887 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 3888 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3889 for (i=0; i<m; i++) { 3890 c->imax[i] = a->imax[i]; 3891 c->ilen[i] = a->ilen[i]; 3892 } 3893 3894 /* allocate the matrix space */ 3895 if (mallocmatspace) { 3896 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 3897 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3898 3899 c->singlemalloc = PETSC_TRUE; 3900 3901 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3902 if (m > 0) { 3903 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3904 if (cpvalues == MAT_COPY_VALUES) { 3905 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3906 } else { 3907 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3908 } 3909 } 3910 } 3911 3912 c->ignorezeroentries = a->ignorezeroentries; 3913 c->roworiented = a->roworiented; 3914 c->nonew = a->nonew; 3915 if (a->diag) { 3916 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 3917 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3918 for (i=0; i<m; i++) { 3919 c->diag[i] = a->diag[i]; 3920 } 3921 } else c->diag = 0; 3922 3923 c->solve_work = 0; 3924 c->saved_values = 0; 3925 c->idiag = 0; 3926 c->ssor_work = 0; 3927 c->keepnonzeropattern = a->keepnonzeropattern; 3928 c->free_a = PETSC_TRUE; 3929 c->free_ij = PETSC_TRUE; 3930 3931 c->rmax = a->rmax; 3932 c->nz = a->nz; 3933 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3934 C->preallocated = PETSC_TRUE; 3935 3936 c->compressedrow.use = a->compressedrow.use; 3937 c->compressedrow.nrows = a->compressedrow.nrows; 3938 if (a->compressedrow.use) { 3939 i = a->compressedrow.nrows; 3940 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 3941 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3942 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3943 } else { 3944 c->compressedrow.use = PETSC_FALSE; 3945 c->compressedrow.i = NULL; 3946 c->compressedrow.rindex = NULL; 3947 } 3948 c->nonzerorowcnt = a->nonzerorowcnt; 3949 C->nonzerostate = A->nonzerostate; 3950 3951 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3952 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3953 PetscFunctionReturn(0); 3954 } 3955 3956 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3957 { 3958 PetscErrorCode ierr; 3959 3960 PetscFunctionBegin; 3961 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 3962 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3963 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 3964 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 3965 } 3966 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 3967 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3968 PetscFunctionReturn(0); 3969 } 3970 3971 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 3972 { 3973 Mat_SeqAIJ *a; 3974 PetscErrorCode ierr; 3975 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3976 int fd; 3977 PetscMPIInt size; 3978 MPI_Comm comm; 3979 PetscInt bs = newMat->rmap->bs; 3980 3981 PetscFunctionBegin; 3982 /* force binary viewer to load .info file if it has not yet done so */ 3983 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 3984 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3985 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3986 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3987 3988 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 3989 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 3990 ierr = PetscOptionsEnd();CHKERRQ(ierr); 3991 if (bs < 0) bs = 1; 3992 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 3993 3994 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3995 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3996 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3997 M = header[1]; N = header[2]; nz = header[3]; 3998 3999 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4000 4001 /* read in row lengths */ 4002 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4003 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4004 4005 /* check if sum of rowlengths is same as nz */ 4006 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4007 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); 4008 4009 /* set global size if not set already*/ 4010 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4011 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4012 } else { 4013 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4014 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4015 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4016 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4017 } 4018 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); 4019 } 4020 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4021 a = (Mat_SeqAIJ*)newMat->data; 4022 4023 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4024 4025 /* read in nonzero values */ 4026 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4027 4028 /* set matrix "i" values */ 4029 a->i[0] = 0; 4030 for (i=1; i<= M; i++) { 4031 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4032 a->ilen[i-1] = rowlengths[i-1]; 4033 } 4034 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4035 4036 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4037 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4038 PetscFunctionReturn(0); 4039 } 4040 4041 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4042 { 4043 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4044 PetscErrorCode ierr; 4045 #if defined(PETSC_USE_COMPLEX) 4046 PetscInt k; 4047 #endif 4048 4049 PetscFunctionBegin; 4050 /* If the matrix dimensions are not equal,or no of nonzeros */ 4051 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4052 *flg = PETSC_FALSE; 4053 PetscFunctionReturn(0); 4054 } 4055 4056 /* if the a->i are the same */ 4057 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4058 if (!*flg) PetscFunctionReturn(0); 4059 4060 /* if a->j are the same */ 4061 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4062 if (!*flg) PetscFunctionReturn(0); 4063 4064 /* if a->a are the same */ 4065 #if defined(PETSC_USE_COMPLEX) 4066 for (k=0; k<a->nz; k++) { 4067 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4068 *flg = PETSC_FALSE; 4069 PetscFunctionReturn(0); 4070 } 4071 } 4072 #else 4073 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4074 #endif 4075 PetscFunctionReturn(0); 4076 } 4077 4078 /*@ 4079 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4080 provided by the user. 4081 4082 Collective on MPI_Comm 4083 4084 Input Parameters: 4085 + comm - must be an MPI communicator of size 1 4086 . m - number of rows 4087 . n - number of columns 4088 . i - row indices 4089 . j - column indices 4090 - a - matrix values 4091 4092 Output Parameter: 4093 . mat - the matrix 4094 4095 Level: intermediate 4096 4097 Notes: 4098 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4099 once the matrix is destroyed and not before 4100 4101 You cannot set new nonzero locations into this matrix, that will generate an error. 4102 4103 The i and j indices are 0 based 4104 4105 The format which is used for the sparse matrix input, is equivalent to a 4106 row-major ordering.. i.e for the following matrix, the input data expected is 4107 as shown 4108 4109 $ 1 0 0 4110 $ 2 0 3 4111 $ 4 5 6 4112 $ 4113 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4114 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4115 $ v = {1,2,3,4,5,6} [size = 6] 4116 4117 4118 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4119 4120 @*/ 4121 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 4122 { 4123 PetscErrorCode ierr; 4124 PetscInt ii; 4125 Mat_SeqAIJ *aij; 4126 #if defined(PETSC_USE_DEBUG) 4127 PetscInt jj; 4128 #endif 4129 4130 PetscFunctionBegin; 4131 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4132 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4133 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4134 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4135 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4136 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4137 aij = (Mat_SeqAIJ*)(*mat)->data; 4138 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4139 4140 aij->i = i; 4141 aij->j = j; 4142 aij->a = a; 4143 aij->singlemalloc = PETSC_FALSE; 4144 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4145 aij->free_a = PETSC_FALSE; 4146 aij->free_ij = PETSC_FALSE; 4147 4148 for (ii=0; ii<m; ii++) { 4149 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4150 #if defined(PETSC_USE_DEBUG) 4151 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]); 4152 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4153 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); 4154 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); 4155 } 4156 #endif 4157 } 4158 #if defined(PETSC_USE_DEBUG) 4159 for (ii=0; ii<aij->i[m]; ii++) { 4160 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4161 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]); 4162 } 4163 #endif 4164 4165 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4166 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4167 PetscFunctionReturn(0); 4168 } 4169 /*@C 4170 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4171 provided by the user. 4172 4173 Collective on MPI_Comm 4174 4175 Input Parameters: 4176 + comm - must be an MPI communicator of size 1 4177 . m - number of rows 4178 . n - number of columns 4179 . i - row indices 4180 . j - column indices 4181 . a - matrix values 4182 . nz - number of nonzeros 4183 - idx - 0 or 1 based 4184 4185 Output Parameter: 4186 . mat - the matrix 4187 4188 Level: intermediate 4189 4190 Notes: 4191 The i and j indices are 0 based 4192 4193 The format which is used for the sparse matrix input, is equivalent to a 4194 row-major ordering.. i.e for the following matrix, the input data expected is 4195 as shown: 4196 4197 1 0 0 4198 2 0 3 4199 4 5 6 4200 4201 i = {0,1,1,2,2,2} 4202 j = {0,0,2,0,1,2} 4203 v = {1,2,3,4,5,6} 4204 4205 4206 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4207 4208 @*/ 4209 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4210 { 4211 PetscErrorCode ierr; 4212 PetscInt ii, *nnz, one = 1,row,col; 4213 4214 4215 PetscFunctionBegin; 4216 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4217 for (ii = 0; ii < nz; ii++) { 4218 nnz[i[ii] - !!idx] += 1; 4219 } 4220 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4221 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4222 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4223 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4224 for (ii = 0; ii < nz; ii++) { 4225 if (idx) { 4226 row = i[ii] - 1; 4227 col = j[ii] - 1; 4228 } else { 4229 row = i[ii]; 4230 col = j[ii]; 4231 } 4232 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4233 } 4234 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4235 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4236 ierr = PetscFree(nnz);CHKERRQ(ierr); 4237 PetscFunctionReturn(0); 4238 } 4239 4240 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4241 { 4242 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4243 PetscErrorCode ierr; 4244 4245 PetscFunctionBegin; 4246 a->idiagvalid = PETSC_FALSE; 4247 a->ibdiagvalid = PETSC_FALSE; 4248 4249 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4250 PetscFunctionReturn(0); 4251 } 4252 4253 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4254 { 4255 PetscErrorCode ierr; 4256 4257 PetscFunctionBegin; 4258 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4259 PetscFunctionReturn(0); 4260 } 4261 4262 /* 4263 Permute A into C's *local* index space using rowemb,colemb. 4264 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4265 of [0,m), colemb is in [0,n). 4266 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4267 */ 4268 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4269 { 4270 /* If making this function public, change the error returned in this function away from _PLIB. */ 4271 PetscErrorCode ierr; 4272 Mat_SeqAIJ *Baij; 4273 PetscBool seqaij; 4274 PetscInt m,n,*nz,i,j,count; 4275 PetscScalar v; 4276 const PetscInt *rowindices,*colindices; 4277 4278 PetscFunctionBegin; 4279 if (!B) PetscFunctionReturn(0); 4280 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4281 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4282 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4283 if (rowemb) { 4284 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4285 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); 4286 } else { 4287 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4288 } 4289 if (colemb) { 4290 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4291 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); 4292 } else { 4293 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4294 } 4295 4296 Baij = (Mat_SeqAIJ*)(B->data); 4297 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4298 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4299 for (i=0; i<B->rmap->n; i++) { 4300 nz[i] = Baij->i[i+1] - Baij->i[i]; 4301 } 4302 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4303 ierr = PetscFree(nz);CHKERRQ(ierr); 4304 } 4305 if (pattern == SUBSET_NONZERO_PATTERN) { 4306 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4307 } 4308 count = 0; 4309 rowindices = NULL; 4310 colindices = NULL; 4311 if (rowemb) { 4312 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4313 } 4314 if (colemb) { 4315 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4316 } 4317 for (i=0; i<B->rmap->n; i++) { 4318 PetscInt row; 4319 row = i; 4320 if (rowindices) row = rowindices[i]; 4321 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4322 PetscInt col; 4323 col = Baij->j[count]; 4324 if (colindices) col = colindices[col]; 4325 v = Baij->a[count]; 4326 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4327 ++count; 4328 } 4329 } 4330 /* FIXME: set C's nonzerostate correctly. */ 4331 /* Assembly for C is necessary. */ 4332 C->preallocated = PETSC_TRUE; 4333 C->assembled = PETSC_TRUE; 4334 C->was_assembled = PETSC_FALSE; 4335 PetscFunctionReturn(0); 4336 } 4337 4338 4339 /* 4340 Special version for direct calls from Fortran 4341 */ 4342 #include <petsc/private/fortranimpl.h> 4343 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4344 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4345 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4346 #define matsetvaluesseqaij_ matsetvaluesseqaij 4347 #endif 4348 4349 /* Change these macros so can be used in void function */ 4350 #undef CHKERRQ 4351 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4352 #undef SETERRQ2 4353 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4354 #undef SETERRQ3 4355 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4356 4357 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) 4358 { 4359 Mat A = *AA; 4360 PetscInt m = *mm, n = *nn; 4361 InsertMode is = *isis; 4362 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4363 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4364 PetscInt *imax,*ai,*ailen; 4365 PetscErrorCode ierr; 4366 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4367 MatScalar *ap,value,*aa; 4368 PetscBool ignorezeroentries = a->ignorezeroentries; 4369 PetscBool roworiented = a->roworiented; 4370 4371 PetscFunctionBegin; 4372 MatCheckPreallocated(A,1); 4373 imax = a->imax; 4374 ai = a->i; 4375 ailen = a->ilen; 4376 aj = a->j; 4377 aa = a->a; 4378 4379 for (k=0; k<m; k++) { /* loop over added rows */ 4380 row = im[k]; 4381 if (row < 0) continue; 4382 #if defined(PETSC_USE_DEBUG) 4383 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4384 #endif 4385 rp = aj + ai[row]; ap = aa + ai[row]; 4386 rmax = imax[row]; nrow = ailen[row]; 4387 low = 0; 4388 high = nrow; 4389 for (l=0; l<n; l++) { /* loop over added columns */ 4390 if (in[l] < 0) continue; 4391 #if defined(PETSC_USE_DEBUG) 4392 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4393 #endif 4394 col = in[l]; 4395 if (roworiented) value = v[l + k*n]; 4396 else value = v[k + l*m]; 4397 4398 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4399 4400 if (col <= lastcol) low = 0; 4401 else high = nrow; 4402 lastcol = col; 4403 while (high-low > 5) { 4404 t = (low+high)/2; 4405 if (rp[t] > col) high = t; 4406 else low = t; 4407 } 4408 for (i=low; i<high; i++) { 4409 if (rp[i] > col) break; 4410 if (rp[i] == col) { 4411 if (is == ADD_VALUES) ap[i] += value; 4412 else ap[i] = value; 4413 goto noinsert; 4414 } 4415 } 4416 if (value == 0.0 && ignorezeroentries) goto noinsert; 4417 if (nonew == 1) goto noinsert; 4418 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4419 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4420 N = nrow++ - 1; a->nz++; high++; 4421 /* shift up all the later entries in this row */ 4422 for (ii=N; ii>=i; ii--) { 4423 rp[ii+1] = rp[ii]; 4424 ap[ii+1] = ap[ii]; 4425 } 4426 rp[i] = col; 4427 ap[i] = value; 4428 A->nonzerostate++; 4429 noinsert:; 4430 low = i + 1; 4431 } 4432 ailen[row] = nrow; 4433 } 4434 PetscFunctionReturnVoid(); 4435 } 4436 4437