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