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