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 /* reset ilen and imax for each row */ 916 for (i=0; i<m; i++) { 917 ailen[i] = imax[i] = ai[i+1] - ai[i]; 918 } 919 a->nz = ai[m]; 920 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); 921 922 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 923 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 924 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 925 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 926 927 A->info.mallocs += a->reallocs; 928 a->reallocs = 0; 929 A->info.nz_unneeded = (double)fshift; 930 a->rmax = rmax; 931 932 ierr = MatCheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 933 934 A->same_nonzero = PETSC_TRUE; 935 936 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 937 938 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 939 PetscFunctionReturn(0); 940 } 941 942 #undef __FUNCT__ 943 #define __FUNCT__ "MatRealPart_SeqAIJ" 944 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 945 { 946 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 947 PetscInt i,nz = a->nz; 948 MatScalar *aa = a->a; 949 PetscErrorCode ierr; 950 951 PetscFunctionBegin; 952 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 953 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 954 PetscFunctionReturn(0); 955 } 956 957 #undef __FUNCT__ 958 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 959 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 960 { 961 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 962 PetscInt i,nz = a->nz; 963 MatScalar *aa = a->a; 964 PetscErrorCode ierr; 965 966 PetscFunctionBegin; 967 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 968 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 969 PetscFunctionReturn(0); 970 } 971 972 #if defined(PETSC_THREADCOMM_ACTIVE) 973 PetscErrorCode MatZeroEntries_SeqAIJ_Kernel(PetscInt thread_id,Mat A) 974 { 975 PetscErrorCode ierr; 976 PetscInt *trstarts=A->rmap->trstarts; 977 PetscInt n,start,end; 978 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 979 980 start = trstarts[thread_id]; 981 end = trstarts[thread_id+1]; 982 n = a->i[end] - a->i[start]; 983 ierr = PetscMemzero(a->a+a->i[start],n*sizeof(PetscScalar));CHKERRQ(ierr); 984 return 0; 985 } 986 987 #undef __FUNCT__ 988 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 989 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 990 { 991 PetscErrorCode ierr; 992 993 PetscFunctionBegin; 994 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatZeroEntries_SeqAIJ_Kernel,1,A);CHKERRQ(ierr); 995 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 996 PetscFunctionReturn(0); 997 } 998 #else 999 #undef __FUNCT__ 1000 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 1001 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1002 { 1003 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1004 PetscErrorCode ierr; 1005 1006 PetscFunctionBegin; 1007 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 1008 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 1009 PetscFunctionReturn(0); 1010 } 1011 #endif 1012 1013 #undef __FUNCT__ 1014 #define __FUNCT__ "MatDestroy_SeqAIJ" 1015 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1016 { 1017 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1018 PetscErrorCode ierr; 1019 1020 PetscFunctionBegin; 1021 #if defined(PETSC_USE_LOG) 1022 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 1023 #endif 1024 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 1025 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 1026 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 1027 ierr = PetscFree(a->diag);CHKERRQ(ierr); 1028 ierr = PetscFree(a->ibdiag);CHKERRQ(ierr); 1029 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 1030 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 1031 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 1032 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 1033 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 1034 ierr = ISColoringDestroy(&a->coloring);CHKERRQ(ierr); 1035 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 1036 ierr = MatDestroy(&a->XtoY);CHKERRQ(ierr); 1037 ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);CHKERRQ(ierr); 1038 ierr = PetscFree(a->matmult_abdense);CHKERRQ(ierr); 1039 1040 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 1041 ierr = PetscFree(A->data);CHKERRQ(ierr); 1042 1043 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1044 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL);CHKERRQ(ierr); 1045 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL);CHKERRQ(ierr); 1046 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL);CHKERRQ(ierr); 1047 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL);CHKERRQ(ierr); 1048 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL);CHKERRQ(ierr); 1049 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL);CHKERRQ(ierr); 1050 ierr = PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL);CHKERRQ(ierr); 1051 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL);CHKERRQ(ierr); 1052 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL);CHKERRQ(ierr); 1053 ierr = PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL);CHKERRQ(ierr); 1054 PetscFunctionReturn(0); 1055 } 1056 1057 #undef __FUNCT__ 1058 #define __FUNCT__ "MatSetOption_SeqAIJ" 1059 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1060 { 1061 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1062 PetscErrorCode ierr; 1063 1064 PetscFunctionBegin; 1065 switch (op) { 1066 case MAT_ROW_ORIENTED: 1067 a->roworiented = flg; 1068 break; 1069 case MAT_KEEP_NONZERO_PATTERN: 1070 a->keepnonzeropattern = flg; 1071 break; 1072 case MAT_NEW_NONZERO_LOCATIONS: 1073 a->nonew = (flg ? 0 : 1); 1074 break; 1075 case MAT_NEW_NONZERO_LOCATION_ERR: 1076 a->nonew = (flg ? -1 : 0); 1077 break; 1078 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1079 a->nonew = (flg ? -2 : 0); 1080 break; 1081 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1082 a->nounused = (flg ? -1 : 0); 1083 break; 1084 case MAT_IGNORE_ZERO_ENTRIES: 1085 a->ignorezeroentries = flg; 1086 break; 1087 case MAT_CHECK_COMPRESSED_ROW: 1088 a->compressedrow.check = 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,nonzerorow=0; 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 nonzerorow += (n>0); 1280 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1281 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1282 y[*ridx++] = sum; 1283 } 1284 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1285 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1286 } else { /* do not use compressed row format */ 1287 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1288 fortranmultaij_(&m,x,ii,aj,aa,y); 1289 #else 1290 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);CHKERRQ(ierr); 1291 #endif 1292 } 1293 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1294 PetscFunctionReturn(0); 1295 } 1296 #else 1297 #undef __FUNCT__ 1298 #define __FUNCT__ "MatMult_SeqAIJ" 1299 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1300 { 1301 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1302 PetscScalar *y; 1303 const PetscScalar *x; 1304 const MatScalar *aa; 1305 PetscErrorCode ierr; 1306 PetscInt m=A->rmap->n; 1307 const PetscInt *aj,*ii,*ridx=NULL; 1308 PetscInt n,i,nonzerorow=0; 1309 PetscScalar sum; 1310 PetscBool usecprow=a->compressedrow.use; 1311 1312 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1313 #pragma disjoint(*x,*y,*aa) 1314 #endif 1315 1316 PetscFunctionBegin; 1317 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1318 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1319 aj = a->j; 1320 aa = a->a; 1321 ii = a->i; 1322 if (usecprow) { /* use compressed row format */ 1323 m = a->compressedrow.nrows; 1324 ii = a->compressedrow.i; 1325 ridx = a->compressedrow.rindex; 1326 for (i=0; i<m; i++) { 1327 n = ii[i+1] - ii[i]; 1328 aj = a->j + ii[i]; 1329 aa = a->a + ii[i]; 1330 sum = 0.0; 1331 nonzerorow += (n>0); 1332 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1333 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1334 y[*ridx++] = sum; 1335 } 1336 } else { /* do not use compressed row format */ 1337 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1338 fortranmultaij_(&m,x,ii,aj,aa,y); 1339 #else 1340 #if defined(PETSC_THREADCOMM_ACTIVE) 1341 ierr = PetscThreadCommRunKernel(PetscObjectComm((PetscObject)A),(PetscThreadKernel)MatMult_SeqAIJ_Kernel,3,A,xx,yy);CHKERRQ(ierr); 1342 #else 1343 for (i=0; i<m; i++) { 1344 n = ii[i+1] - ii[i]; 1345 aj = a->j + ii[i]; 1346 aa = a->a + ii[i]; 1347 sum = 0.0; 1348 nonzerorow += (n>0); 1349 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1350 y[i] = sum; 1351 } 1352 #endif 1353 #endif 1354 } 1355 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1356 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1357 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1358 PetscFunctionReturn(0); 1359 } 1360 #endif 1361 1362 #undef __FUNCT__ 1363 #define __FUNCT__ "MatMultMax_SeqAIJ" 1364 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1365 { 1366 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1367 PetscScalar *y; 1368 const PetscScalar *x; 1369 const MatScalar *aa; 1370 PetscErrorCode ierr; 1371 PetscInt m=A->rmap->n; 1372 const PetscInt *aj,*ii,*ridx=NULL; 1373 PetscInt n,i,nonzerorow=0; 1374 PetscScalar sum; 1375 PetscBool usecprow=a->compressedrow.use; 1376 1377 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1378 #pragma disjoint(*x,*y,*aa) 1379 #endif 1380 1381 PetscFunctionBegin; 1382 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1383 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1384 aj = a->j; 1385 aa = a->a; 1386 ii = a->i; 1387 if (usecprow) { /* use compressed row format */ 1388 m = a->compressedrow.nrows; 1389 ii = a->compressedrow.i; 1390 ridx = a->compressedrow.rindex; 1391 for (i=0; i<m; i++) { 1392 n = ii[i+1] - ii[i]; 1393 aj = a->j + ii[i]; 1394 aa = a->a + ii[i]; 1395 sum = 0.0; 1396 nonzerorow += (n>0); 1397 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1398 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1399 y[*ridx++] = sum; 1400 } 1401 } else { /* do not use compressed row format */ 1402 for (i=0; i<m; i++) { 1403 n = ii[i+1] - ii[i]; 1404 aj = a->j + ii[i]; 1405 aa = a->a + ii[i]; 1406 sum = 0.0; 1407 nonzerorow += (n>0); 1408 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1409 y[i] = sum; 1410 } 1411 } 1412 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1413 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1414 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1415 PetscFunctionReturn(0); 1416 } 1417 1418 #undef __FUNCT__ 1419 #define __FUNCT__ "MatMultAddMax_SeqAIJ" 1420 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1421 { 1422 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1423 PetscScalar *y,*z; 1424 const PetscScalar *x; 1425 const MatScalar *aa; 1426 PetscErrorCode ierr; 1427 PetscInt m = A->rmap->n,*aj,*ii; 1428 PetscInt n,i,*ridx=NULL; 1429 PetscScalar sum; 1430 PetscBool usecprow=a->compressedrow.use; 1431 1432 PetscFunctionBegin; 1433 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1434 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1435 if (zz != yy) { 1436 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1437 } else { 1438 z = y; 1439 } 1440 1441 aj = a->j; 1442 aa = a->a; 1443 ii = a->i; 1444 if (usecprow) { /* use compressed row format */ 1445 if (zz != yy) { 1446 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1447 } 1448 m = a->compressedrow.nrows; 1449 ii = a->compressedrow.i; 1450 ridx = a->compressedrow.rindex; 1451 for (i=0; i<m; i++) { 1452 n = ii[i+1] - ii[i]; 1453 aj = a->j + ii[i]; 1454 aa = a->a + ii[i]; 1455 sum = y[*ridx]; 1456 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1457 z[*ridx++] = sum; 1458 } 1459 } else { /* do not use compressed row format */ 1460 for (i=0; i<m; i++) { 1461 n = ii[i+1] - ii[i]; 1462 aj = a->j + ii[i]; 1463 aa = a->a + ii[i]; 1464 sum = y[i]; 1465 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1466 z[i] = sum; 1467 } 1468 } 1469 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1470 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1471 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1472 if (zz != yy) { 1473 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1474 } 1475 PetscFunctionReturn(0); 1476 } 1477 1478 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1479 #undef __FUNCT__ 1480 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1481 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1482 { 1483 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1484 PetscScalar *y,*z; 1485 const PetscScalar *x; 1486 const MatScalar *aa; 1487 PetscErrorCode ierr; 1488 PetscInt m = A->rmap->n,*aj,*ii; 1489 PetscInt n,i,*ridx=NULL; 1490 PetscScalar sum; 1491 PetscBool usecprow=a->compressedrow.use; 1492 1493 PetscFunctionBegin; 1494 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1495 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1496 if (zz != yy) { 1497 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1498 } else { 1499 z = y; 1500 } 1501 1502 aj = a->j; 1503 aa = a->a; 1504 ii = a->i; 1505 if (usecprow) { /* use compressed row format */ 1506 if (zz != yy) { 1507 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1508 } 1509 m = a->compressedrow.nrows; 1510 ii = a->compressedrow.i; 1511 ridx = a->compressedrow.rindex; 1512 for (i=0; i<m; i++) { 1513 n = ii[i+1] - ii[i]; 1514 aj = a->j + ii[i]; 1515 aa = a->a + ii[i]; 1516 sum = y[*ridx]; 1517 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1518 z[*ridx++] = sum; 1519 } 1520 } else { /* do not use compressed row format */ 1521 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1522 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1523 #else 1524 for (i=0; i<m; i++) { 1525 n = ii[i+1] - ii[i]; 1526 aj = a->j + ii[i]; 1527 aa = a->a + ii[i]; 1528 sum = y[i]; 1529 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1530 z[i] = sum; 1531 } 1532 #endif 1533 } 1534 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1535 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1536 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1537 if (zz != yy) { 1538 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1539 } 1540 #if defined(PETSC_HAVE_CUSP) 1541 /* 1542 ierr = VecView(xx,0);CHKERRQ(ierr); 1543 ierr = VecView(zz,0);CHKERRQ(ierr); 1544 ierr = MatView(A,0);CHKERRQ(ierr); 1545 */ 1546 #endif 1547 PetscFunctionReturn(0); 1548 } 1549 1550 /* 1551 Adds diagonal pointers to sparse matrix structure. 1552 */ 1553 #undef __FUNCT__ 1554 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1555 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1556 { 1557 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1558 PetscErrorCode ierr; 1559 PetscInt i,j,m = A->rmap->n; 1560 1561 PetscFunctionBegin; 1562 if (!a->diag) { 1563 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1564 ierr = PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));CHKERRQ(ierr); 1565 } 1566 for (i=0; i<A->rmap->n; i++) { 1567 a->diag[i] = a->i[i+1]; 1568 for (j=a->i[i]; j<a->i[i+1]; j++) { 1569 if (a->j[j] == i) { 1570 a->diag[i] = j; 1571 break; 1572 } 1573 } 1574 } 1575 PetscFunctionReturn(0); 1576 } 1577 1578 /* 1579 Checks for missing diagonals 1580 */ 1581 #undef __FUNCT__ 1582 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1583 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1584 { 1585 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1586 PetscInt *diag,*jj = a->j,i; 1587 1588 PetscFunctionBegin; 1589 *missing = PETSC_FALSE; 1590 if (A->rmap->n > 0 && !jj) { 1591 *missing = PETSC_TRUE; 1592 if (d) *d = 0; 1593 PetscInfo(A,"Matrix has no entries therefore is missing diagonal"); 1594 } else { 1595 diag = a->diag; 1596 for (i=0; i<A->rmap->n; i++) { 1597 if (jj[diag[i]] != i) { 1598 *missing = PETSC_TRUE; 1599 if (d) *d = i; 1600 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1601 break; 1602 } 1603 } 1604 } 1605 PetscFunctionReturn(0); 1606 } 1607 1608 #undef __FUNCT__ 1609 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1610 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1611 { 1612 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1613 PetscErrorCode ierr; 1614 PetscInt i,*diag,m = A->rmap->n; 1615 MatScalar *v = a->a; 1616 PetscScalar *idiag,*mdiag; 1617 1618 PetscFunctionBegin; 1619 if (a->idiagvalid) PetscFunctionReturn(0); 1620 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1621 diag = a->diag; 1622 if (!a->idiag) { 1623 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1624 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1625 v = a->a; 1626 } 1627 mdiag = a->mdiag; 1628 idiag = a->idiag; 1629 1630 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1631 for (i=0; i<m; i++) { 1632 mdiag[i] = v[diag[i]]; 1633 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1634 idiag[i] = 1.0/v[diag[i]]; 1635 } 1636 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1637 } else { 1638 for (i=0; i<m; i++) { 1639 mdiag[i] = v[diag[i]]; 1640 idiag[i] = omega/(fshift + v[diag[i]]); 1641 } 1642 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1643 } 1644 a->idiagvalid = PETSC_TRUE; 1645 PetscFunctionReturn(0); 1646 } 1647 1648 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1649 #undef __FUNCT__ 1650 #define __FUNCT__ "MatSOR_SeqAIJ" 1651 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1652 { 1653 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1654 PetscScalar *x,d,sum,*t,scale; 1655 const MatScalar *v = a->a,*idiag=0,*mdiag; 1656 const PetscScalar *b, *bs,*xb, *ts; 1657 PetscErrorCode ierr; 1658 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1659 const PetscInt *idx,*diag; 1660 1661 PetscFunctionBegin; 1662 its = its*lits; 1663 1664 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1665 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1666 a->fshift = fshift; 1667 a->omega = omega; 1668 1669 diag = a->diag; 1670 t = a->ssor_work; 1671 idiag = a->idiag; 1672 mdiag = a->mdiag; 1673 1674 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1675 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1676 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1677 if (flag == SOR_APPLY_UPPER) { 1678 /* apply (U + D/omega) to the vector */ 1679 bs = b; 1680 for (i=0; i<m; i++) { 1681 d = fshift + mdiag[i]; 1682 n = a->i[i+1] - diag[i] - 1; 1683 idx = a->j + diag[i] + 1; 1684 v = a->a + diag[i] + 1; 1685 sum = b[i]*d/omega; 1686 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1687 x[i] = sum; 1688 } 1689 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1690 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1691 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1692 PetscFunctionReturn(0); 1693 } 1694 1695 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1696 else if (flag & SOR_EISENSTAT) { 1697 /* Let A = L + U + D; where L is lower trianglar, 1698 U is upper triangular, E = D/omega; This routine applies 1699 1700 (L + E)^{-1} A (U + E)^{-1} 1701 1702 to a vector efficiently using Eisenstat's trick. 1703 */ 1704 scale = (2.0/omega) - 1.0; 1705 1706 /* x = (E + U)^{-1} b */ 1707 for (i=m-1; i>=0; i--) { 1708 n = a->i[i+1] - diag[i] - 1; 1709 idx = a->j + diag[i] + 1; 1710 v = a->a + diag[i] + 1; 1711 sum = b[i]; 1712 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1713 x[i] = sum*idiag[i]; 1714 } 1715 1716 /* t = b - (2*E - D)x */ 1717 v = a->a; 1718 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1719 1720 /* t = (E + L)^{-1}t */ 1721 ts = t; 1722 diag = a->diag; 1723 for (i=0; i<m; i++) { 1724 n = diag[i] - a->i[i]; 1725 idx = a->j + a->i[i]; 1726 v = a->a + a->i[i]; 1727 sum = t[i]; 1728 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1729 t[i] = sum*idiag[i]; 1730 /* x = x + t */ 1731 x[i] += t[i]; 1732 } 1733 1734 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1735 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1736 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1737 PetscFunctionReturn(0); 1738 } 1739 if (flag & SOR_ZERO_INITIAL_GUESS) { 1740 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1741 for (i=0; i<m; i++) { 1742 n = diag[i] - a->i[i]; 1743 idx = a->j + a->i[i]; 1744 v = a->a + a->i[i]; 1745 sum = b[i]; 1746 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1747 t[i] = sum; 1748 x[i] = sum*idiag[i]; 1749 } 1750 xb = t; 1751 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1752 } else xb = b; 1753 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1754 for (i=m-1; i>=0; i--) { 1755 n = a->i[i+1] - diag[i] - 1; 1756 idx = a->j + diag[i] + 1; 1757 v = a->a + diag[i] + 1; 1758 sum = xb[i]; 1759 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1760 if (xb == b) { 1761 x[i] = sum*idiag[i]; 1762 } else { 1763 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1764 } 1765 } 1766 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1767 } 1768 its--; 1769 } 1770 while (its--) { 1771 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1772 for (i=0; i<m; i++) { 1773 /* lower */ 1774 n = diag[i] - a->i[i]; 1775 idx = a->j + a->i[i]; 1776 v = a->a + a->i[i]; 1777 sum = b[i]; 1778 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1779 t[i] = sum; /* save application of the lower-triangular part */ 1780 /* upper */ 1781 n = a->i[i+1] - diag[i] - 1; 1782 idx = a->j + diag[i] + 1; 1783 v = a->a + diag[i] + 1; 1784 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1785 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1786 } 1787 xb = t; 1788 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1789 } else xb = b; 1790 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1791 for (i=m-1; i>=0; i--) { 1792 sum = xb[i]; 1793 if (xb == b) { 1794 /* whole matrix (no checkpointing available) */ 1795 n = a->i[i+1] - a->i[i]; 1796 idx = a->j + a->i[i]; 1797 v = a->a + a->i[i]; 1798 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1799 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1800 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1801 n = a->i[i+1] - diag[i] - 1; 1802 idx = a->j + diag[i] + 1; 1803 v = a->a + diag[i] + 1; 1804 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1805 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1806 } 1807 } 1808 if (xb == b) { 1809 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1810 } else { 1811 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1812 } 1813 } 1814 } 1815 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1816 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1817 PetscFunctionReturn(0); 1818 } 1819 1820 1821 #undef __FUNCT__ 1822 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1823 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1824 { 1825 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1826 1827 PetscFunctionBegin; 1828 info->block_size = 1.0; 1829 info->nz_allocated = (double)a->maxnz; 1830 info->nz_used = (double)a->nz; 1831 info->nz_unneeded = (double)(a->maxnz - a->nz); 1832 info->assemblies = (double)A->num_ass; 1833 info->mallocs = (double)A->info.mallocs; 1834 info->memory = ((PetscObject)A)->mem; 1835 if (A->factortype) { 1836 info->fill_ratio_given = A->info.fill_ratio_given; 1837 info->fill_ratio_needed = A->info.fill_ratio_needed; 1838 info->factor_mallocs = A->info.factor_mallocs; 1839 } else { 1840 info->fill_ratio_given = 0; 1841 info->fill_ratio_needed = 0; 1842 info->factor_mallocs = 0; 1843 } 1844 PetscFunctionReturn(0); 1845 } 1846 1847 #undef __FUNCT__ 1848 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1849 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1850 { 1851 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1852 PetscInt i,m = A->rmap->n - 1,d = 0; 1853 PetscErrorCode ierr; 1854 const PetscScalar *xx; 1855 PetscScalar *bb; 1856 PetscBool missing; 1857 1858 PetscFunctionBegin; 1859 if (x && b) { 1860 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1861 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1862 for (i=0; i<N; i++) { 1863 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1864 bb[rows[i]] = diag*xx[rows[i]]; 1865 } 1866 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1867 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1868 } 1869 1870 if (a->keepnonzeropattern) { 1871 for (i=0; i<N; i++) { 1872 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1873 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1874 } 1875 if (diag != 0.0) { 1876 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1877 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1878 for (i=0; i<N; i++) { 1879 a->a[a->diag[rows[i]]] = diag; 1880 } 1881 } 1882 A->same_nonzero = PETSC_TRUE; 1883 } else { 1884 if (diag != 0.0) { 1885 for (i=0; i<N; i++) { 1886 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1887 if (a->ilen[rows[i]] > 0) { 1888 a->ilen[rows[i]] = 1; 1889 a->a[a->i[rows[i]]] = diag; 1890 a->j[a->i[rows[i]]] = rows[i]; 1891 } else { /* in case row was completely empty */ 1892 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1893 } 1894 } 1895 } else { 1896 for (i=0; i<N; i++) { 1897 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1898 a->ilen[rows[i]] = 0; 1899 } 1900 } 1901 A->same_nonzero = PETSC_FALSE; 1902 } 1903 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1904 PetscFunctionReturn(0); 1905 } 1906 1907 #undef __FUNCT__ 1908 #define __FUNCT__ "MatZeroRowsColumns_SeqAIJ" 1909 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1910 { 1911 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1912 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1913 PetscErrorCode ierr; 1914 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 1915 const PetscScalar *xx; 1916 PetscScalar *bb; 1917 1918 PetscFunctionBegin; 1919 if (x && b) { 1920 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1921 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1922 vecs = PETSC_TRUE; 1923 } 1924 ierr = PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);CHKERRQ(ierr); 1925 ierr = PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));CHKERRQ(ierr); 1926 for (i=0; i<N; i++) { 1927 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1928 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1929 1930 zeroed[rows[i]] = PETSC_TRUE; 1931 } 1932 for (i=0; i<A->rmap->n; i++) { 1933 if (!zeroed[i]) { 1934 for (j=a->i[i]; j<a->i[i+1]; j++) { 1935 if (zeroed[a->j[j]]) { 1936 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 1937 a->a[j] = 0.0; 1938 } 1939 } 1940 } else if (vecs) bb[i] = diag*xx[i]; 1941 } 1942 if (x && b) { 1943 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1944 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1945 } 1946 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1947 if (diag != 0.0) { 1948 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1949 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1950 for (i=0; i<N; i++) { 1951 a->a[a->diag[rows[i]]] = diag; 1952 } 1953 } 1954 A->same_nonzero = PETSC_TRUE; 1955 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1956 PetscFunctionReturn(0); 1957 } 1958 1959 #undef __FUNCT__ 1960 #define __FUNCT__ "MatGetRow_SeqAIJ" 1961 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1962 { 1963 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1964 PetscInt *itmp; 1965 1966 PetscFunctionBegin; 1967 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1968 1969 *nz = a->i[row+1] - a->i[row]; 1970 if (v) *v = a->a + a->i[row]; 1971 if (idx) { 1972 itmp = a->j + a->i[row]; 1973 if (*nz) *idx = itmp; 1974 else *idx = 0; 1975 } 1976 PetscFunctionReturn(0); 1977 } 1978 1979 /* remove this function? */ 1980 #undef __FUNCT__ 1981 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1982 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1983 { 1984 PetscFunctionBegin; 1985 PetscFunctionReturn(0); 1986 } 1987 1988 #undef __FUNCT__ 1989 #define __FUNCT__ "MatNorm_SeqAIJ" 1990 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1991 { 1992 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1993 MatScalar *v = a->a; 1994 PetscReal sum = 0.0; 1995 PetscErrorCode ierr; 1996 PetscInt i,j; 1997 1998 PetscFunctionBegin; 1999 if (type == NORM_FROBENIUS) { 2000 for (i=0; i<a->nz; i++) { 2001 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2002 } 2003 *nrm = PetscSqrtReal(sum); 2004 } else if (type == NORM_1) { 2005 PetscReal *tmp; 2006 PetscInt *jj = a->j; 2007 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 2008 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 2009 *nrm = 0.0; 2010 for (j=0; j<a->nz; j++) { 2011 tmp[*jj++] += PetscAbsScalar(*v); v++; 2012 } 2013 for (j=0; j<A->cmap->n; j++) { 2014 if (tmp[j] > *nrm) *nrm = tmp[j]; 2015 } 2016 ierr = PetscFree(tmp);CHKERRQ(ierr); 2017 } else if (type == NORM_INFINITY) { 2018 *nrm = 0.0; 2019 for (j=0; j<A->rmap->n; j++) { 2020 v = a->a + a->i[j]; 2021 sum = 0.0; 2022 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2023 sum += PetscAbsScalar(*v); v++; 2024 } 2025 if (sum > *nrm) *nrm = sum; 2026 } 2027 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2028 PetscFunctionReturn(0); 2029 } 2030 2031 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2032 #undef __FUNCT__ 2033 #define __FUNCT__ "MatTransposeSymbolic_SeqAIJ" 2034 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2035 { 2036 PetscErrorCode ierr; 2037 PetscInt i,j,anzj; 2038 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2039 PetscInt an=A->cmap->N,am=A->rmap->N; 2040 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2041 2042 PetscFunctionBegin; 2043 /* Allocate space for symbolic transpose info and work array */ 2044 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 2045 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 2046 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 2047 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 2048 2049 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2050 /* Note: offset by 1 for fast conversion into csr format. */ 2051 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2052 /* Form ati for csr format of A^T. */ 2053 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2054 2055 /* Copy ati into atfill so we have locations of the next free space in atj */ 2056 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2057 2058 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2059 for (i=0;i<am;i++) { 2060 anzj = ai[i+1] - ai[i]; 2061 for (j=0;j<anzj;j++) { 2062 atj[atfill[*aj]] = i; 2063 atfill[*aj++] += 1; 2064 } 2065 } 2066 2067 /* Clean up temporary space and complete requests. */ 2068 ierr = PetscFree(atfill);CHKERRQ(ierr); 2069 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2070 2071 (*B)->rmap->bs = A->cmap->bs; 2072 (*B)->cmap->bs = A->rmap->bs; 2073 2074 b = (Mat_SeqAIJ*)((*B)->data); 2075 b->free_a = PETSC_FALSE; 2076 b->free_ij = PETSC_TRUE; 2077 b->nonew = 0; 2078 PetscFunctionReturn(0); 2079 } 2080 2081 #undef __FUNCT__ 2082 #define __FUNCT__ "MatTranspose_SeqAIJ" 2083 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2084 { 2085 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2086 Mat C; 2087 PetscErrorCode ierr; 2088 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2089 MatScalar *array = a->a; 2090 2091 PetscFunctionBegin; 2092 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"); 2093 2094 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 2095 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 2096 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 2097 2098 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2099 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2100 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2101 ierr = MatSetBlockSizes(C,A->cmap->bs,A->rmap->bs);CHKERRQ(ierr); 2102 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2103 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2104 ierr = PetscFree(col);CHKERRQ(ierr); 2105 } else { 2106 C = *B; 2107 } 2108 2109 for (i=0; i<m; i++) { 2110 len = ai[i+1]-ai[i]; 2111 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2112 array += len; 2113 aj += len; 2114 } 2115 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2116 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2117 2118 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 2119 *B = C; 2120 } else { 2121 ierr = MatHeaderMerge(A,C);CHKERRQ(ierr); 2122 } 2123 PetscFunctionReturn(0); 2124 } 2125 2126 #undef __FUNCT__ 2127 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 2128 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2129 { 2130 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2131 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2132 MatScalar *va,*vb; 2133 PetscErrorCode ierr; 2134 PetscInt ma,na,mb,nb, i; 2135 2136 PetscFunctionBegin; 2137 bij = (Mat_SeqAIJ*) B->data; 2138 2139 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2140 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2141 if (ma!=nb || na!=mb) { 2142 *f = PETSC_FALSE; 2143 PetscFunctionReturn(0); 2144 } 2145 aii = aij->i; bii = bij->i; 2146 adx = aij->j; bdx = bij->j; 2147 va = aij->a; vb = bij->a; 2148 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 2149 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 2150 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2151 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2152 2153 *f = PETSC_TRUE; 2154 for (i=0; i<ma; i++) { 2155 while (aptr[i]<aii[i+1]) { 2156 PetscInt idc,idr; 2157 PetscScalar vc,vr; 2158 /* column/row index/value */ 2159 idc = adx[aptr[i]]; 2160 idr = bdx[bptr[idc]]; 2161 vc = va[aptr[i]]; 2162 vr = vb[bptr[idc]]; 2163 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2164 *f = PETSC_FALSE; 2165 goto done; 2166 } else { 2167 aptr[i]++; 2168 if (B || i!=idc) bptr[idc]++; 2169 } 2170 } 2171 } 2172 done: 2173 ierr = PetscFree(aptr);CHKERRQ(ierr); 2174 ierr = PetscFree(bptr);CHKERRQ(ierr); 2175 PetscFunctionReturn(0); 2176 } 2177 2178 #undef __FUNCT__ 2179 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 2180 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2181 { 2182 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) A->data; 2183 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2184 MatScalar *va,*vb; 2185 PetscErrorCode ierr; 2186 PetscInt ma,na,mb,nb, i; 2187 2188 PetscFunctionBegin; 2189 bij = (Mat_SeqAIJ*) B->data; 2190 2191 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2192 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2193 if (ma!=nb || na!=mb) { 2194 *f = PETSC_FALSE; 2195 PetscFunctionReturn(0); 2196 } 2197 aii = aij->i; bii = bij->i; 2198 adx = aij->j; bdx = bij->j; 2199 va = aij->a; vb = bij->a; 2200 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 2201 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 2202 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2203 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2204 2205 *f = PETSC_TRUE; 2206 for (i=0; i<ma; i++) { 2207 while (aptr[i]<aii[i+1]) { 2208 PetscInt idc,idr; 2209 PetscScalar vc,vr; 2210 /* column/row index/value */ 2211 idc = adx[aptr[i]]; 2212 idr = bdx[bptr[idc]]; 2213 vc = va[aptr[i]]; 2214 vr = vb[bptr[idc]]; 2215 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2216 *f = PETSC_FALSE; 2217 goto done; 2218 } else { 2219 aptr[i]++; 2220 if (B || i!=idc) bptr[idc]++; 2221 } 2222 } 2223 } 2224 done: 2225 ierr = PetscFree(aptr);CHKERRQ(ierr); 2226 ierr = PetscFree(bptr);CHKERRQ(ierr); 2227 PetscFunctionReturn(0); 2228 } 2229 2230 #undef __FUNCT__ 2231 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 2232 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2233 { 2234 PetscErrorCode ierr; 2235 2236 PetscFunctionBegin; 2237 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2238 PetscFunctionReturn(0); 2239 } 2240 2241 #undef __FUNCT__ 2242 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 2243 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2244 { 2245 PetscErrorCode ierr; 2246 2247 PetscFunctionBegin; 2248 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2249 PetscFunctionReturn(0); 2250 } 2251 2252 #undef __FUNCT__ 2253 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 2254 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2255 { 2256 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2257 PetscScalar *l,*r,x; 2258 MatScalar *v; 2259 PetscErrorCode ierr; 2260 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 2261 2262 PetscFunctionBegin; 2263 if (ll) { 2264 /* The local size is used so that VecMPI can be passed to this routine 2265 by MatDiagonalScale_MPIAIJ */ 2266 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2267 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2268 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 2269 v = a->a; 2270 for (i=0; i<m; i++) { 2271 x = l[i]; 2272 M = a->i[i+1] - a->i[i]; 2273 for (j=0; j<M; j++) (*v++) *= x; 2274 } 2275 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 2276 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2277 } 2278 if (rr) { 2279 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2280 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2281 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 2282 v = a->a; jj = a->j; 2283 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2284 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 2285 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2286 } 2287 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2288 PetscFunctionReturn(0); 2289 } 2290 2291 #undef __FUNCT__ 2292 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 2293 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2294 { 2295 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2296 PetscErrorCode ierr; 2297 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2298 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2299 const PetscInt *irow,*icol; 2300 PetscInt nrows,ncols; 2301 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2302 MatScalar *a_new,*mat_a; 2303 Mat C; 2304 PetscBool stride,sorted; 2305 2306 PetscFunctionBegin; 2307 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 2308 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 2309 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 2310 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 2311 2312 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2313 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2314 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2315 2316 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2317 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2318 if (stride && step == 1) { 2319 /* special case of contiguous rows */ 2320 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 2321 /* loop over new rows determining lens and starting points */ 2322 for (i=0; i<nrows; i++) { 2323 kstart = ai[irow[i]]; 2324 kend = kstart + ailen[irow[i]]; 2325 for (k=kstart; k<kend; k++) { 2326 if (aj[k] >= first) { 2327 starts[i] = k; 2328 break; 2329 } 2330 } 2331 sum = 0; 2332 while (k < kend) { 2333 if (aj[k++] >= first+ncols) break; 2334 sum++; 2335 } 2336 lens[i] = sum; 2337 } 2338 /* create submatrix */ 2339 if (scall == MAT_REUSE_MATRIX) { 2340 PetscInt n_cols,n_rows; 2341 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2342 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2343 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2344 C = *B; 2345 } else { 2346 PetscInt rbs,cbs; 2347 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2348 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2349 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2350 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2351 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2352 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2353 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2354 } 2355 c = (Mat_SeqAIJ*)C->data; 2356 2357 /* loop over rows inserting into submatrix */ 2358 a_new = c->a; 2359 j_new = c->j; 2360 i_new = c->i; 2361 2362 for (i=0; i<nrows; i++) { 2363 ii = starts[i]; 2364 lensi = lens[i]; 2365 for (k=0; k<lensi; k++) { 2366 *j_new++ = aj[ii+k] - first; 2367 } 2368 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2369 a_new += lensi; 2370 i_new[i+1] = i_new[i] + lensi; 2371 c->ilen[i] = lensi; 2372 } 2373 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2374 } else { 2375 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2376 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 2377 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 2378 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2379 for (i=0; i<ncols; i++) { 2380 #if defined(PETSC_USE_DEBUG) 2381 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); 2382 #endif 2383 smap[icol[i]] = i+1; 2384 } 2385 2386 /* determine lens of each row */ 2387 for (i=0; i<nrows; i++) { 2388 kstart = ai[irow[i]]; 2389 kend = kstart + a->ilen[irow[i]]; 2390 lens[i] = 0; 2391 for (k=kstart; k<kend; k++) { 2392 if (smap[aj[k]]) { 2393 lens[i]++; 2394 } 2395 } 2396 } 2397 /* Create and fill new matrix */ 2398 if (scall == MAT_REUSE_MATRIX) { 2399 PetscBool equal; 2400 2401 c = (Mat_SeqAIJ*)((*B)->data); 2402 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2403 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2404 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2405 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2406 C = *B; 2407 } else { 2408 PetscInt rbs,cbs; 2409 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2410 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2411 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2412 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2413 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2414 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2415 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2416 } 2417 c = (Mat_SeqAIJ*)(C->data); 2418 for (i=0; i<nrows; i++) { 2419 row = irow[i]; 2420 kstart = ai[row]; 2421 kend = kstart + a->ilen[row]; 2422 mat_i = c->i[i]; 2423 mat_j = c->j + mat_i; 2424 mat_a = c->a + mat_i; 2425 mat_ilen = c->ilen + i; 2426 for (k=kstart; k<kend; k++) { 2427 if ((tcol=smap[a->j[k]])) { 2428 *mat_j++ = tcol - 1; 2429 *mat_a++ = a->a[k]; 2430 (*mat_ilen)++; 2431 2432 } 2433 } 2434 } 2435 /* Free work space */ 2436 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2437 ierr = PetscFree(smap);CHKERRQ(ierr); 2438 ierr = PetscFree(lens);CHKERRQ(ierr); 2439 } 2440 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2441 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2442 2443 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2444 *B = C; 2445 PetscFunctionReturn(0); 2446 } 2447 2448 #undef __FUNCT__ 2449 #define __FUNCT__ "MatGetMultiProcBlock_SeqAIJ" 2450 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2451 { 2452 PetscErrorCode ierr; 2453 Mat B; 2454 2455 PetscFunctionBegin; 2456 if (scall == MAT_INITIAL_MATRIX) { 2457 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2458 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2459 ierr = MatSetBlockSizes(B,mat->rmap->bs,mat->cmap->bs);CHKERRQ(ierr); 2460 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2461 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2462 *subMat = B; 2463 } else { 2464 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2465 } 2466 PetscFunctionReturn(0); 2467 } 2468 2469 #undef __FUNCT__ 2470 #define __FUNCT__ "MatILUFactor_SeqAIJ" 2471 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2472 { 2473 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2474 PetscErrorCode ierr; 2475 Mat outA; 2476 PetscBool row_identity,col_identity; 2477 2478 PetscFunctionBegin; 2479 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2480 2481 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2482 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2483 2484 outA = inA; 2485 outA->factortype = MAT_FACTOR_LU; 2486 2487 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2488 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2489 2490 a->row = row; 2491 2492 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2493 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2494 2495 a->col = col; 2496 2497 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2498 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2499 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2500 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2501 2502 if (!a->solve_work) { /* this matrix may have been factored before */ 2503 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 2504 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2505 } 2506 2507 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2508 if (row_identity && col_identity) { 2509 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2510 } else { 2511 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2512 } 2513 PetscFunctionReturn(0); 2514 } 2515 2516 #undef __FUNCT__ 2517 #define __FUNCT__ "MatScale_SeqAIJ" 2518 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2519 { 2520 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2521 PetscScalar oalpha = alpha; 2522 PetscErrorCode ierr; 2523 PetscBLASInt one = 1,bnz; 2524 2525 PetscFunctionBegin; 2526 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2527 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2528 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2529 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2530 PetscFunctionReturn(0); 2531 } 2532 2533 #undef __FUNCT__ 2534 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 2535 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2536 { 2537 PetscErrorCode ierr; 2538 PetscInt i; 2539 2540 PetscFunctionBegin; 2541 if (scall == MAT_INITIAL_MATRIX) { 2542 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 2543 } 2544 2545 for (i=0; i<n; i++) { 2546 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2547 } 2548 PetscFunctionReturn(0); 2549 } 2550 2551 #undef __FUNCT__ 2552 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 2553 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2554 { 2555 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2556 PetscErrorCode ierr; 2557 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2558 const PetscInt *idx; 2559 PetscInt start,end,*ai,*aj; 2560 PetscBT table; 2561 2562 PetscFunctionBegin; 2563 m = A->rmap->n; 2564 ai = a->i; 2565 aj = a->j; 2566 2567 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2568 2569 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 2570 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2571 2572 for (i=0; i<is_max; i++) { 2573 /* Initialize the two local arrays */ 2574 isz = 0; 2575 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2576 2577 /* Extract the indices, assume there can be duplicate entries */ 2578 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2579 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2580 2581 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2582 for (j=0; j<n; ++j) { 2583 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2584 } 2585 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2586 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2587 2588 k = 0; 2589 for (j=0; j<ov; j++) { /* for each overlap */ 2590 n = isz; 2591 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2592 row = nidx[k]; 2593 start = ai[row]; 2594 end = ai[row+1]; 2595 for (l = start; l<end; l++) { 2596 val = aj[l]; 2597 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2598 } 2599 } 2600 } 2601 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2602 } 2603 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2604 ierr = PetscFree(nidx);CHKERRQ(ierr); 2605 PetscFunctionReturn(0); 2606 } 2607 2608 /* -------------------------------------------------------------- */ 2609 #undef __FUNCT__ 2610 #define __FUNCT__ "MatPermute_SeqAIJ" 2611 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2612 { 2613 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2614 PetscErrorCode ierr; 2615 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2616 const PetscInt *row,*col; 2617 PetscInt *cnew,j,*lens; 2618 IS icolp,irowp; 2619 PetscInt *cwork = NULL; 2620 PetscScalar *vwork = NULL; 2621 2622 PetscFunctionBegin; 2623 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2624 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2625 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2626 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2627 2628 /* determine lengths of permuted rows */ 2629 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2630 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2631 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2632 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2633 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 2634 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2635 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2636 ierr = PetscFree(lens);CHKERRQ(ierr); 2637 2638 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2639 for (i=0; i<m; i++) { 2640 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2641 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2642 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2643 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2644 } 2645 ierr = PetscFree(cnew);CHKERRQ(ierr); 2646 2647 (*B)->assembled = PETSC_FALSE; 2648 2649 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2650 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2651 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2652 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2653 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2654 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2655 PetscFunctionReturn(0); 2656 } 2657 2658 #undef __FUNCT__ 2659 #define __FUNCT__ "MatCopy_SeqAIJ" 2660 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2661 { 2662 PetscErrorCode ierr; 2663 2664 PetscFunctionBegin; 2665 /* If the two matrices have the same copy implementation, use fast copy. */ 2666 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2667 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2668 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2669 2670 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"); 2671 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2672 } else { 2673 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2674 } 2675 PetscFunctionReturn(0); 2676 } 2677 2678 #undef __FUNCT__ 2679 #define __FUNCT__ "MatSetUp_SeqAIJ" 2680 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2681 { 2682 PetscErrorCode ierr; 2683 2684 PetscFunctionBegin; 2685 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2686 PetscFunctionReturn(0); 2687 } 2688 2689 #undef __FUNCT__ 2690 #define __FUNCT__ "MatSeqAIJGetArray_SeqAIJ" 2691 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2692 { 2693 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2694 2695 PetscFunctionBegin; 2696 *array = a->a; 2697 PetscFunctionReturn(0); 2698 } 2699 2700 #undef __FUNCT__ 2701 #define __FUNCT__ "MatSeqAIJRestoreArray_SeqAIJ" 2702 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2703 { 2704 PetscFunctionBegin; 2705 PetscFunctionReturn(0); 2706 } 2707 2708 /* 2709 Computes the number of nonzeros per row needed for preallocation when X and Y 2710 have different nonzero structure. 2711 */ 2712 #undef __FUNCT__ 2713 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2714 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2715 { 2716 PetscInt i,m=Y->rmap->N; 2717 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2718 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2719 const PetscInt *xi = x->i,*yi = y->i; 2720 2721 PetscFunctionBegin; 2722 /* Set the number of nonzeros in the new matrix */ 2723 for (i=0; i<m; i++) { 2724 PetscInt j,k,nzx = xi[i+1] - xi[i],nzy = yi[i+1] - yi[i]; 2725 const PetscInt *xj = x->j+xi[i],*yj = y->j+yi[i]; 2726 nnz[i] = 0; 2727 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2728 for (; k<nzy && yj[k]<xj[j]; k++) nnz[i]++; /* Catch up to X */ 2729 if (k<nzy && yj[k]==xj[j]) k++; /* Skip duplicate */ 2730 nnz[i]++; 2731 } 2732 for (; k<nzy; k++) nnz[i]++; 2733 } 2734 PetscFunctionReturn(0); 2735 } 2736 2737 #undef __FUNCT__ 2738 #define __FUNCT__ "MatAXPY_SeqAIJ" 2739 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2740 { 2741 PetscErrorCode ierr; 2742 PetscInt i; 2743 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2744 PetscBLASInt one=1,bnz; 2745 2746 PetscFunctionBegin; 2747 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2748 if (str == SAME_NONZERO_PATTERN) { 2749 PetscScalar alpha = a; 2750 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2751 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2752 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2753 if (y->xtoy && y->XtoY != X) { 2754 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2755 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 2756 } 2757 if (!y->xtoy) { /* get xtoy */ 2758 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,NULL, y->i,y->j,NULL, &y->xtoy);CHKERRQ(ierr); 2759 y->XtoY = X; 2760 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2761 } 2762 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2763 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); 2764 } else { 2765 Mat B; 2766 PetscInt *nnz; 2767 ierr = PetscMalloc(Y->rmap->N*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 2768 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2769 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2770 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2771 ierr = MatSetBlockSizes(B,Y->rmap->bs,Y->cmap->bs);CHKERRQ(ierr); 2772 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2773 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2774 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2775 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2776 ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr); 2777 ierr = PetscFree(nnz);CHKERRQ(ierr); 2778 } 2779 PetscFunctionReturn(0); 2780 } 2781 2782 #undef __FUNCT__ 2783 #define __FUNCT__ "MatConjugate_SeqAIJ" 2784 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2785 { 2786 #if defined(PETSC_USE_COMPLEX) 2787 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2788 PetscInt i,nz; 2789 PetscScalar *a; 2790 2791 PetscFunctionBegin; 2792 nz = aij->nz; 2793 a = aij->a; 2794 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2795 #else 2796 PetscFunctionBegin; 2797 #endif 2798 PetscFunctionReturn(0); 2799 } 2800 2801 #undef __FUNCT__ 2802 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2803 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2804 { 2805 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2806 PetscErrorCode ierr; 2807 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2808 PetscReal atmp; 2809 PetscScalar *x; 2810 MatScalar *aa; 2811 2812 PetscFunctionBegin; 2813 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2814 aa = a->a; 2815 ai = a->i; 2816 aj = a->j; 2817 2818 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2819 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2820 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2821 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2822 for (i=0; i<m; i++) { 2823 ncols = ai[1] - ai[0]; ai++; 2824 x[i] = 0.0; 2825 for (j=0; j<ncols; j++) { 2826 atmp = PetscAbsScalar(*aa); 2827 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2828 aa++; aj++; 2829 } 2830 } 2831 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2832 PetscFunctionReturn(0); 2833 } 2834 2835 #undef __FUNCT__ 2836 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2837 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2838 { 2839 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2840 PetscErrorCode ierr; 2841 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2842 PetscScalar *x; 2843 MatScalar *aa; 2844 2845 PetscFunctionBegin; 2846 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2847 aa = a->a; 2848 ai = a->i; 2849 aj = a->j; 2850 2851 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2852 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2853 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2854 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2855 for (i=0; i<m; i++) { 2856 ncols = ai[1] - ai[0]; ai++; 2857 if (ncols == A->cmap->n) { /* row is dense */ 2858 x[i] = *aa; if (idx) idx[i] = 0; 2859 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2860 x[i] = 0.0; 2861 if (idx) { 2862 idx[i] = 0; /* in case ncols is zero */ 2863 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2864 if (aj[j] > j) { 2865 idx[i] = j; 2866 break; 2867 } 2868 } 2869 } 2870 } 2871 for (j=0; j<ncols; j++) { 2872 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2873 aa++; aj++; 2874 } 2875 } 2876 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2877 PetscFunctionReturn(0); 2878 } 2879 2880 #undef __FUNCT__ 2881 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2882 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2883 { 2884 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2885 PetscErrorCode ierr; 2886 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2887 PetscReal atmp; 2888 PetscScalar *x; 2889 MatScalar *aa; 2890 2891 PetscFunctionBegin; 2892 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2893 aa = a->a; 2894 ai = a->i; 2895 aj = a->j; 2896 2897 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2898 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2899 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2900 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); 2901 for (i=0; i<m; i++) { 2902 ncols = ai[1] - ai[0]; ai++; 2903 if (ncols) { 2904 /* Get first nonzero */ 2905 for (j = 0; j < ncols; j++) { 2906 atmp = PetscAbsScalar(aa[j]); 2907 if (atmp > 1.0e-12) { 2908 x[i] = atmp; 2909 if (idx) idx[i] = aj[j]; 2910 break; 2911 } 2912 } 2913 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2914 } else { 2915 x[i] = 0.0; if (idx) idx[i] = 0; 2916 } 2917 for (j = 0; j < ncols; j++) { 2918 atmp = PetscAbsScalar(*aa); 2919 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2920 aa++; aj++; 2921 } 2922 } 2923 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2924 PetscFunctionReturn(0); 2925 } 2926 2927 #undef __FUNCT__ 2928 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2929 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2930 { 2931 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2932 PetscErrorCode ierr; 2933 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2934 PetscScalar *x; 2935 MatScalar *aa; 2936 2937 PetscFunctionBegin; 2938 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2939 aa = a->a; 2940 ai = a->i; 2941 aj = a->j; 2942 2943 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2944 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2945 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2946 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2947 for (i=0; i<m; i++) { 2948 ncols = ai[1] - ai[0]; ai++; 2949 if (ncols == A->cmap->n) { /* row is dense */ 2950 x[i] = *aa; if (idx) idx[i] = 0; 2951 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2952 x[i] = 0.0; 2953 if (idx) { /* find first implicit 0.0 in the row */ 2954 idx[i] = 0; /* in case ncols is zero */ 2955 for (j=0; j<ncols; j++) { 2956 if (aj[j] > j) { 2957 idx[i] = j; 2958 break; 2959 } 2960 } 2961 } 2962 } 2963 for (j=0; j<ncols; j++) { 2964 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2965 aa++; aj++; 2966 } 2967 } 2968 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2969 PetscFunctionReturn(0); 2970 } 2971 2972 #include <petscblaslapack.h> 2973 #include <petsc-private/kernels/blockinvert.h> 2974 2975 #undef __FUNCT__ 2976 #define __FUNCT__ "MatInvertBlockDiagonal_SeqAIJ" 2977 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 2978 { 2979 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 2980 PetscErrorCode ierr; 2981 PetscInt i,bs = A->rmap->bs,mbs = A->rmap->n/A->rmap->bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 2982 MatScalar *diag,work[25],*v_work; 2983 PetscReal shift = 0.0; 2984 2985 PetscFunctionBegin; 2986 if (a->ibdiagvalid) { 2987 if (values) *values = a->ibdiag; 2988 PetscFunctionReturn(0); 2989 } 2990 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 2991 if (!a->ibdiag) { 2992 ierr = PetscMalloc(bs2*mbs*sizeof(PetscScalar),&a->ibdiag);CHKERRQ(ierr); 2993 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2994 } 2995 diag = a->ibdiag; 2996 if (values) *values = a->ibdiag; 2997 /* factor and invert each block */ 2998 switch (bs) { 2999 case 1: 3000 for (i=0; i<mbs; i++) { 3001 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3002 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3003 } 3004 break; 3005 case 2: 3006 for (i=0; i<mbs; i++) { 3007 ij[0] = 2*i; ij[1] = 2*i + 1; 3008 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3009 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift);CHKERRQ(ierr); 3010 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3011 diag += 4; 3012 } 3013 break; 3014 case 3: 3015 for (i=0; i<mbs; i++) { 3016 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3017 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3018 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift);CHKERRQ(ierr); 3019 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3020 diag += 9; 3021 } 3022 break; 3023 case 4: 3024 for (i=0; i<mbs; i++) { 3025 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3026 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3027 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift);CHKERRQ(ierr); 3028 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3029 diag += 16; 3030 } 3031 break; 3032 case 5: 3033 for (i=0; i<mbs; i++) { 3034 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3035 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3036 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift);CHKERRQ(ierr); 3037 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3038 diag += 25; 3039 } 3040 break; 3041 case 6: 3042 for (i=0; i<mbs; i++) { 3043 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; 3044 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3045 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift);CHKERRQ(ierr); 3046 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3047 diag += 36; 3048 } 3049 break; 3050 case 7: 3051 for (i=0; i<mbs; i++) { 3052 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; 3053 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3054 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift);CHKERRQ(ierr); 3055 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3056 diag += 49; 3057 } 3058 break; 3059 default: 3060 ierr = PetscMalloc3(bs,MatScalar,&v_work,bs,PetscInt,&v_pivots,bs,PetscInt,&IJ);CHKERRQ(ierr); 3061 for (i=0; i<mbs; i++) { 3062 for (j=0; j<bs; j++) { 3063 IJ[j] = bs*i + j; 3064 } 3065 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3066 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);CHKERRQ(ierr); 3067 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3068 diag += bs2; 3069 } 3070 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3071 } 3072 a->ibdiagvalid = PETSC_TRUE; 3073 PetscFunctionReturn(0); 3074 } 3075 3076 #undef __FUNCT__ 3077 #define __FUNCT__ "MatSetRandom_SeqAIJ" 3078 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3079 { 3080 PetscErrorCode ierr; 3081 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3082 PetscScalar a; 3083 PetscInt m,n,i,j,col; 3084 3085 PetscFunctionBegin; 3086 if (!x->assembled) { 3087 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3088 for (i=0; i<m; i++) { 3089 for (j=0; j<aij->imax[i]; j++) { 3090 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3091 col = (PetscInt)(n*PetscRealPart(a)); 3092 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3093 } 3094 } 3095 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3096 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3097 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3098 PetscFunctionReturn(0); 3099 } 3100 3101 /* -------------------------------------------------------------------*/ 3102 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3103 MatGetRow_SeqAIJ, 3104 MatRestoreRow_SeqAIJ, 3105 MatMult_SeqAIJ, 3106 /* 4*/ MatMultAdd_SeqAIJ, 3107 MatMultTranspose_SeqAIJ, 3108 MatMultTransposeAdd_SeqAIJ, 3109 0, 3110 0, 3111 0, 3112 /* 10*/ 0, 3113 MatLUFactor_SeqAIJ, 3114 0, 3115 MatSOR_SeqAIJ, 3116 MatTranspose_SeqAIJ, 3117 /*1 5*/ MatGetInfo_SeqAIJ, 3118 MatEqual_SeqAIJ, 3119 MatGetDiagonal_SeqAIJ, 3120 MatDiagonalScale_SeqAIJ, 3121 MatNorm_SeqAIJ, 3122 /* 20*/ 0, 3123 MatAssemblyEnd_SeqAIJ, 3124 MatSetOption_SeqAIJ, 3125 MatZeroEntries_SeqAIJ, 3126 /* 24*/ MatZeroRows_SeqAIJ, 3127 0, 3128 0, 3129 0, 3130 0, 3131 /* 29*/ MatSetUp_SeqAIJ, 3132 0, 3133 0, 3134 0, 3135 0, 3136 /* 34*/ MatDuplicate_SeqAIJ, 3137 0, 3138 0, 3139 MatILUFactor_SeqAIJ, 3140 0, 3141 /* 39*/ MatAXPY_SeqAIJ, 3142 MatGetSubMatrices_SeqAIJ, 3143 MatIncreaseOverlap_SeqAIJ, 3144 MatGetValues_SeqAIJ, 3145 MatCopy_SeqAIJ, 3146 /* 44*/ MatGetRowMax_SeqAIJ, 3147 MatScale_SeqAIJ, 3148 0, 3149 MatDiagonalSet_SeqAIJ, 3150 MatZeroRowsColumns_SeqAIJ, 3151 /* 49*/ MatSetRandom_SeqAIJ, 3152 MatGetRowIJ_SeqAIJ, 3153 MatRestoreRowIJ_SeqAIJ, 3154 MatGetColumnIJ_SeqAIJ, 3155 MatRestoreColumnIJ_SeqAIJ, 3156 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3157 0, 3158 0, 3159 MatPermute_SeqAIJ, 3160 0, 3161 /* 59*/ 0, 3162 MatDestroy_SeqAIJ, 3163 MatView_SeqAIJ, 3164 0, 3165 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3166 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3167 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3168 0, 3169 0, 3170 0, 3171 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3172 MatGetRowMinAbs_SeqAIJ, 3173 0, 3174 MatSetColoring_SeqAIJ, 3175 0, 3176 /* 74*/ MatSetValuesAdifor_SeqAIJ, 3177 MatFDColoringApply_AIJ, 3178 0, 3179 0, 3180 0, 3181 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3182 0, 3183 0, 3184 0, 3185 MatLoad_SeqAIJ, 3186 /* 84*/ MatIsSymmetric_SeqAIJ, 3187 MatIsHermitian_SeqAIJ, 3188 0, 3189 0, 3190 0, 3191 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3192 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3193 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3194 MatPtAP_SeqAIJ_SeqAIJ, 3195 MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy, 3196 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3197 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3198 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3199 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3200 0, 3201 /* 99*/ 0, 3202 0, 3203 0, 3204 MatConjugate_SeqAIJ, 3205 0, 3206 /*104*/ MatSetValuesRow_SeqAIJ, 3207 MatRealPart_SeqAIJ, 3208 MatImaginaryPart_SeqAIJ, 3209 0, 3210 0, 3211 /*109*/ MatMatSolve_SeqAIJ, 3212 0, 3213 MatGetRowMin_SeqAIJ, 3214 0, 3215 MatMissingDiagonal_SeqAIJ, 3216 /*114*/ 0, 3217 0, 3218 0, 3219 0, 3220 0, 3221 /*119*/ 0, 3222 0, 3223 0, 3224 0, 3225 MatGetMultiProcBlock_SeqAIJ, 3226 /*124*/ MatFindNonzeroRows_SeqAIJ, 3227 MatGetColumnNorms_SeqAIJ, 3228 MatInvertBlockDiagonal_SeqAIJ, 3229 0, 3230 0, 3231 /*129*/ 0, 3232 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3233 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3234 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3235 MatTransposeColoringCreate_SeqAIJ, 3236 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3237 MatTransColoringApplyDenToSp_SeqAIJ, 3238 MatRARt_SeqAIJ_SeqAIJ, 3239 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3240 MatRARtNumeric_SeqAIJ_SeqAIJ, 3241 /*139*/0, 3242 0 3243 }; 3244 3245 #undef __FUNCT__ 3246 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3247 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3248 { 3249 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3250 PetscInt i,nz,n; 3251 3252 PetscFunctionBegin; 3253 nz = aij->maxnz; 3254 n = mat->rmap->n; 3255 for (i=0; i<nz; i++) { 3256 aij->j[i] = indices[i]; 3257 } 3258 aij->nz = nz; 3259 for (i=0; i<n; i++) { 3260 aij->ilen[i] = aij->imax[i]; 3261 } 3262 PetscFunctionReturn(0); 3263 } 3264 3265 #undef __FUNCT__ 3266 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3267 /*@ 3268 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3269 in the matrix. 3270 3271 Input Parameters: 3272 + mat - the SeqAIJ matrix 3273 - indices - the column indices 3274 3275 Level: advanced 3276 3277 Notes: 3278 This can be called if you have precomputed the nonzero structure of the 3279 matrix and want to provide it to the matrix object to improve the performance 3280 of the MatSetValues() operation. 3281 3282 You MUST have set the correct numbers of nonzeros per row in the call to 3283 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3284 3285 MUST be called before any calls to MatSetValues(); 3286 3287 The indices should start with zero, not one. 3288 3289 @*/ 3290 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3291 { 3292 PetscErrorCode ierr; 3293 3294 PetscFunctionBegin; 3295 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3296 PetscValidPointer(indices,2); 3297 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3298 PetscFunctionReturn(0); 3299 } 3300 3301 /* ----------------------------------------------------------------------------------------*/ 3302 3303 #undef __FUNCT__ 3304 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3305 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3306 { 3307 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3308 PetscErrorCode ierr; 3309 size_t nz = aij->i[mat->rmap->n]; 3310 3311 PetscFunctionBegin; 3312 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3313 3314 /* allocate space for values if not already there */ 3315 if (!aij->saved_values) { 3316 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 3317 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3318 } 3319 3320 /* copy values over */ 3321 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3322 PetscFunctionReturn(0); 3323 } 3324 3325 #undef __FUNCT__ 3326 #define __FUNCT__ "MatStoreValues" 3327 /*@ 3328 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3329 example, reuse of the linear part of a Jacobian, while recomputing the 3330 nonlinear portion. 3331 3332 Collect on Mat 3333 3334 Input Parameters: 3335 . mat - the matrix (currently only AIJ matrices support this option) 3336 3337 Level: advanced 3338 3339 Common Usage, with SNESSolve(): 3340 $ Create Jacobian matrix 3341 $ Set linear terms into matrix 3342 $ Apply boundary conditions to matrix, at this time matrix must have 3343 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3344 $ boundary conditions again will not change the nonzero structure 3345 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3346 $ ierr = MatStoreValues(mat); 3347 $ Call SNESSetJacobian() with matrix 3348 $ In your Jacobian routine 3349 $ ierr = MatRetrieveValues(mat); 3350 $ Set nonlinear terms in matrix 3351 3352 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3353 $ // build linear portion of Jacobian 3354 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3355 $ ierr = MatStoreValues(mat); 3356 $ loop over nonlinear iterations 3357 $ ierr = MatRetrieveValues(mat); 3358 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3359 $ // call MatAssemblyBegin/End() on matrix 3360 $ Solve linear system with Jacobian 3361 $ endloop 3362 3363 Notes: 3364 Matrix must already be assemblied before calling this routine 3365 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3366 calling this routine. 3367 3368 When this is called multiple times it overwrites the previous set of stored values 3369 and does not allocated additional space. 3370 3371 .seealso: MatRetrieveValues() 3372 3373 @*/ 3374 PetscErrorCode MatStoreValues(Mat mat) 3375 { 3376 PetscErrorCode ierr; 3377 3378 PetscFunctionBegin; 3379 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3380 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3381 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3382 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3383 PetscFunctionReturn(0); 3384 } 3385 3386 #undef __FUNCT__ 3387 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3388 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3389 { 3390 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3391 PetscErrorCode ierr; 3392 PetscInt nz = aij->i[mat->rmap->n]; 3393 3394 PetscFunctionBegin; 3395 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3396 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3397 /* copy values over */ 3398 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3399 PetscFunctionReturn(0); 3400 } 3401 3402 #undef __FUNCT__ 3403 #define __FUNCT__ "MatRetrieveValues" 3404 /*@ 3405 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3406 example, reuse of the linear part of a Jacobian, while recomputing the 3407 nonlinear portion. 3408 3409 Collect on Mat 3410 3411 Input Parameters: 3412 . mat - the matrix (currently on AIJ matrices support this option) 3413 3414 Level: advanced 3415 3416 .seealso: MatStoreValues() 3417 3418 @*/ 3419 PetscErrorCode MatRetrieveValues(Mat mat) 3420 { 3421 PetscErrorCode ierr; 3422 3423 PetscFunctionBegin; 3424 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3425 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3426 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3427 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3428 PetscFunctionReturn(0); 3429 } 3430 3431 3432 /* --------------------------------------------------------------------------------*/ 3433 #undef __FUNCT__ 3434 #define __FUNCT__ "MatCreateSeqAIJ" 3435 /*@C 3436 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3437 (the default parallel PETSc format). For good matrix assembly performance 3438 the user should preallocate the matrix storage by setting the parameter nz 3439 (or the array nnz). By setting these parameters accurately, performance 3440 during matrix assembly can be increased by more than a factor of 50. 3441 3442 Collective on MPI_Comm 3443 3444 Input Parameters: 3445 + comm - MPI communicator, set to PETSC_COMM_SELF 3446 . m - number of rows 3447 . n - number of columns 3448 . nz - number of nonzeros per row (same for all rows) 3449 - nnz - array containing the number of nonzeros in the various rows 3450 (possibly different for each row) or NULL 3451 3452 Output Parameter: 3453 . A - the matrix 3454 3455 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3456 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3457 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3458 3459 Notes: 3460 If nnz is given then nz is ignored 3461 3462 The AIJ format (also called the Yale sparse matrix format or 3463 compressed row storage), is fully compatible with standard Fortran 77 3464 storage. That is, the stored row and column indices can begin at 3465 either one (as in Fortran) or zero. See the users' manual for details. 3466 3467 Specify the preallocated storage with either nz or nnz (not both). 3468 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3469 allocation. For large problems you MUST preallocate memory or you 3470 will get TERRIBLE performance, see the users' manual chapter on matrices. 3471 3472 By default, this format uses inodes (identical nodes) when possible, to 3473 improve numerical efficiency of matrix-vector products and solves. We 3474 search for consecutive rows with the same nonzero structure, thereby 3475 reusing matrix information to achieve increased efficiency. 3476 3477 Options Database Keys: 3478 + -mat_no_inode - Do not use inodes 3479 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3480 3481 Level: intermediate 3482 3483 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3484 3485 @*/ 3486 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3487 { 3488 PetscErrorCode ierr; 3489 3490 PetscFunctionBegin; 3491 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3492 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3493 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3494 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3495 PetscFunctionReturn(0); 3496 } 3497 3498 #undef __FUNCT__ 3499 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3500 /*@C 3501 MatSeqAIJSetPreallocation - For good matrix assembly performance 3502 the user should preallocate the matrix storage by setting the parameter nz 3503 (or the array nnz). By setting these parameters accurately, performance 3504 during matrix assembly can be increased by more than a factor of 50. 3505 3506 Collective on MPI_Comm 3507 3508 Input Parameters: 3509 + B - The matrix-free 3510 . nz - number of nonzeros per row (same for all rows) 3511 - nnz - array containing the number of nonzeros in the various rows 3512 (possibly different for each row) or NULL 3513 3514 Notes: 3515 If nnz is given then nz is ignored 3516 3517 The AIJ format (also called the Yale sparse matrix format or 3518 compressed row storage), is fully compatible with standard Fortran 77 3519 storage. That is, the stored row and column indices can begin at 3520 either one (as in Fortran) or zero. See the users' manual for details. 3521 3522 Specify the preallocated storage with either nz or nnz (not both). 3523 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3524 allocation. For large problems you MUST preallocate memory or you 3525 will get TERRIBLE performance, see the users' manual chapter on matrices. 3526 3527 You can call MatGetInfo() to get information on how effective the preallocation was; 3528 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3529 You can also run with the option -info and look for messages with the string 3530 malloc in them to see if additional memory allocation was needed. 3531 3532 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3533 entries or columns indices 3534 3535 By default, this format uses inodes (identical nodes) when possible, to 3536 improve numerical efficiency of matrix-vector products and solves. We 3537 search for consecutive rows with the same nonzero structure, thereby 3538 reusing matrix information to achieve increased efficiency. 3539 3540 Options Database Keys: 3541 + -mat_no_inode - Do not use inodes 3542 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3543 - -mat_aij_oneindex - Internally use indexing starting at 1 3544 rather than 0. Note that when calling MatSetValues(), 3545 the user still MUST index entries starting at 0! 3546 3547 Level: intermediate 3548 3549 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3550 3551 @*/ 3552 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3553 { 3554 PetscErrorCode ierr; 3555 3556 PetscFunctionBegin; 3557 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3558 PetscValidType(B,1); 3559 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3560 PetscFunctionReturn(0); 3561 } 3562 3563 #undef __FUNCT__ 3564 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3565 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3566 { 3567 Mat_SeqAIJ *b; 3568 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3569 PetscErrorCode ierr; 3570 PetscInt i; 3571 3572 PetscFunctionBegin; 3573 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3574 if (nz == MAT_SKIP_ALLOCATION) { 3575 skipallocation = PETSC_TRUE; 3576 nz = 0; 3577 } 3578 3579 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3580 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3581 3582 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3583 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3584 if (nnz) { 3585 for (i=0; i<B->rmap->n; i++) { 3586 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]); 3587 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); 3588 } 3589 } 3590 3591 B->preallocated = PETSC_TRUE; 3592 3593 b = (Mat_SeqAIJ*)B->data; 3594 3595 if (!skipallocation) { 3596 if (!b->imax) { 3597 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3598 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3599 } 3600 if (!nnz) { 3601 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3602 else if (nz < 0) nz = 1; 3603 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3604 nz = nz*B->rmap->n; 3605 } else { 3606 nz = 0; 3607 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3608 } 3609 /* b->ilen will count nonzeros in each row so far. */ 3610 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3611 3612 /* allocate the matrix space */ 3613 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3614 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3615 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3616 b->i[0] = 0; 3617 for (i=1; i<B->rmap->n+1; i++) { 3618 b->i[i] = b->i[i-1] + b->imax[i-1]; 3619 } 3620 b->singlemalloc = PETSC_TRUE; 3621 b->free_a = PETSC_TRUE; 3622 b->free_ij = PETSC_TRUE; 3623 #if defined(PETSC_THREADCOMM_ACTIVE) 3624 ierr = MatZeroEntries_SeqAIJ(B);CHKERRQ(ierr); 3625 #endif 3626 } else { 3627 b->free_a = PETSC_FALSE; 3628 b->free_ij = PETSC_FALSE; 3629 } 3630 3631 b->nz = 0; 3632 b->maxnz = nz; 3633 B->info.nz_unneeded = (double)b->maxnz; 3634 if (realalloc) { 3635 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3636 } 3637 PetscFunctionReturn(0); 3638 } 3639 3640 #undef __FUNCT__ 3641 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3642 /*@ 3643 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3644 3645 Input Parameters: 3646 + B - the matrix 3647 . i - the indices into j for the start of each row (starts with zero) 3648 . j - the column indices for each row (starts with zero) these must be sorted for each row 3649 - v - optional values in the matrix 3650 3651 Level: developer 3652 3653 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3654 3655 .keywords: matrix, aij, compressed row, sparse, sequential 3656 3657 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3658 @*/ 3659 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3660 { 3661 PetscErrorCode ierr; 3662 3663 PetscFunctionBegin; 3664 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3665 PetscValidType(B,1); 3666 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3667 PetscFunctionReturn(0); 3668 } 3669 3670 #undef __FUNCT__ 3671 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3672 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3673 { 3674 PetscInt i; 3675 PetscInt m,n; 3676 PetscInt nz; 3677 PetscInt *nnz, nz_max = 0; 3678 PetscScalar *values; 3679 PetscErrorCode ierr; 3680 3681 PetscFunctionBegin; 3682 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3683 3684 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3685 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3686 3687 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3688 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3689 for (i = 0; i < m; i++) { 3690 nz = Ii[i+1]- Ii[i]; 3691 nz_max = PetscMax(nz_max, nz); 3692 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3693 nnz[i] = nz; 3694 } 3695 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3696 ierr = PetscFree(nnz);CHKERRQ(ierr); 3697 3698 if (v) { 3699 values = (PetscScalar*) v; 3700 } else { 3701 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3702 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3703 } 3704 3705 for (i = 0; i < m; i++) { 3706 nz = Ii[i+1] - Ii[i]; 3707 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3708 } 3709 3710 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3711 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3712 3713 if (!v) { 3714 ierr = PetscFree(values);CHKERRQ(ierr); 3715 } 3716 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3717 PetscFunctionReturn(0); 3718 } 3719 3720 #include <../src/mat/impls/dense/seq/dense.h> 3721 #include <petsc-private/kernels/petscaxpy.h> 3722 3723 #undef __FUNCT__ 3724 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3725 /* 3726 Computes (B'*A')' since computing B*A directly is untenable 3727 3728 n p p 3729 ( ) ( ) ( ) 3730 m ( A ) * n ( B ) = m ( C ) 3731 ( ) ( ) ( ) 3732 3733 */ 3734 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3735 { 3736 PetscErrorCode ierr; 3737 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3738 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3739 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3740 PetscInt i,n,m,q,p; 3741 const PetscInt *ii,*idx; 3742 const PetscScalar *b,*a,*a_q; 3743 PetscScalar *c,*c_q; 3744 3745 PetscFunctionBegin; 3746 m = A->rmap->n; 3747 n = A->cmap->n; 3748 p = B->cmap->n; 3749 a = sub_a->v; 3750 b = sub_b->a; 3751 c = sub_c->v; 3752 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3753 3754 ii = sub_b->i; 3755 idx = sub_b->j; 3756 for (i=0; i<n; i++) { 3757 q = ii[i+1] - ii[i]; 3758 while (q-->0) { 3759 c_q = c + m*(*idx); 3760 a_q = a + m*i; 3761 PetscKernelAXPY(c_q,*b,a_q,m); 3762 idx++; 3763 b++; 3764 } 3765 } 3766 PetscFunctionReturn(0); 3767 } 3768 3769 #undef __FUNCT__ 3770 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3771 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3772 { 3773 PetscErrorCode ierr; 3774 PetscInt m=A->rmap->n,n=B->cmap->n; 3775 Mat Cmat; 3776 3777 PetscFunctionBegin; 3778 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); 3779 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3780 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3781 ierr = MatSetBlockSizes(Cmat,A->rmap->bs,B->cmap->bs);CHKERRQ(ierr); 3782 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3783 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3784 3785 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3786 3787 *C = Cmat; 3788 PetscFunctionReturn(0); 3789 } 3790 3791 /* ----------------------------------------------------------------*/ 3792 #undef __FUNCT__ 3793 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3794 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3795 { 3796 PetscErrorCode ierr; 3797 3798 PetscFunctionBegin; 3799 if (scall == MAT_INITIAL_MATRIX) { 3800 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3801 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3802 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3803 } 3804 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3805 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3806 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3807 PetscFunctionReturn(0); 3808 } 3809 3810 3811 /*MC 3812 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3813 based on compressed sparse row format. 3814 3815 Options Database Keys: 3816 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3817 3818 Level: beginner 3819 3820 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3821 M*/ 3822 3823 /*MC 3824 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3825 3826 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3827 and MATMPIAIJ otherwise. As a result, for single process communicators, 3828 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3829 for communicators controlling multiple processes. It is recommended that you call both of 3830 the above preallocation routines for simplicity. 3831 3832 Options Database Keys: 3833 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3834 3835 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3836 enough exist. 3837 3838 Level: beginner 3839 3840 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3841 M*/ 3842 3843 /*MC 3844 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3845 3846 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3847 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3848 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3849 for communicators controlling multiple processes. It is recommended that you call both of 3850 the above preallocation routines for simplicity. 3851 3852 Options Database Keys: 3853 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3854 3855 Level: beginner 3856 3857 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3858 M*/ 3859 3860 #if defined(PETSC_HAVE_PASTIX) 3861 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3862 #endif 3863 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3864 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat*); 3865 #endif 3866 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3867 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3868 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3869 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool*); 3870 #if defined(PETSC_HAVE_MUMPS) 3871 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3872 #endif 3873 #if defined(PETSC_HAVE_SUPERLU) 3874 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3875 #endif 3876 #if defined(PETSC_HAVE_SUPERLU_DIST) 3877 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3878 #endif 3879 #if defined(PETSC_HAVE_UMFPACK) 3880 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3881 #endif 3882 #if defined(PETSC_HAVE_CHOLMOD) 3883 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3884 #endif 3885 #if defined(PETSC_HAVE_LUSOL) 3886 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3887 #endif 3888 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3889 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3890 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3891 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3892 #endif 3893 #if defined(PETSC_HAVE_CLIQUE) 3894 PETSC_EXTERN PetscErrorCode MatGetFactor_aij_clique(Mat,MatFactorType,Mat*); 3895 #endif 3896 3897 3898 #undef __FUNCT__ 3899 #define __FUNCT__ "MatSeqAIJGetArray" 3900 /*@C 3901 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 3902 3903 Not Collective 3904 3905 Input Parameter: 3906 . mat - a MATSEQDENSE matrix 3907 3908 Output Parameter: 3909 . array - pointer to the data 3910 3911 Level: intermediate 3912 3913 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3914 @*/ 3915 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3916 { 3917 PetscErrorCode ierr; 3918 3919 PetscFunctionBegin; 3920 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3921 PetscFunctionReturn(0); 3922 } 3923 3924 #undef __FUNCT__ 3925 #define __FUNCT__ "MatSeqAIJRestoreArray" 3926 /*@C 3927 MatSeqAIJRestoreArray - returns access to the array where the data for a SeqSeqAIJ matrix is stored obtained by MatSeqAIJGetArray() 3928 3929 Not Collective 3930 3931 Input Parameters: 3932 . mat - a MATSEQDENSE matrix 3933 . array - pointer to the data 3934 3935 Level: intermediate 3936 3937 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 3938 @*/ 3939 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3940 { 3941 PetscErrorCode ierr; 3942 3943 PetscFunctionBegin; 3944 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3945 PetscFunctionReturn(0); 3946 } 3947 3948 #undef __FUNCT__ 3949 #define __FUNCT__ "MatCreate_SeqAIJ" 3950 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 3951 { 3952 Mat_SeqAIJ *b; 3953 PetscErrorCode ierr; 3954 PetscMPIInt size; 3955 3956 PetscFunctionBegin; 3957 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3958 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3959 3960 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3961 3962 B->data = (void*)b; 3963 3964 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3965 3966 b->row = 0; 3967 b->col = 0; 3968 b->icol = 0; 3969 b->reallocs = 0; 3970 b->ignorezeroentries = PETSC_FALSE; 3971 b->roworiented = PETSC_TRUE; 3972 b->nonew = 0; 3973 b->diag = 0; 3974 b->solve_work = 0; 3975 B->spptr = 0; 3976 b->saved_values = 0; 3977 b->idiag = 0; 3978 b->mdiag = 0; 3979 b->ssor_work = 0; 3980 b->omega = 1.0; 3981 b->fshift = 0.0; 3982 b->idiagvalid = PETSC_FALSE; 3983 b->ibdiagvalid = PETSC_FALSE; 3984 b->keepnonzeropattern = PETSC_FALSE; 3985 b->xtoy = 0; 3986 b->XtoY = 0; 3987 B->same_nonzero = PETSC_FALSE; 3988 3989 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3990 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 3991 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 3992 3993 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3994 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_matlab_C",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3995 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3996 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3997 #endif 3998 #if defined(PETSC_HAVE_PASTIX) 3999 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_pastix_C",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 4000 #endif 4001 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 4002 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_essl_C",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 4003 #endif 4004 #if defined(PETSC_HAVE_SUPERLU) 4005 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_C",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 4006 #endif 4007 #if defined(PETSC_HAVE_SUPERLU_DIST) 4008 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_superlu_dist_C",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 4009 #endif 4010 #if defined(PETSC_HAVE_MUMPS) 4011 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_mumps_C",MatGetFactor_aij_mumps);CHKERRQ(ierr); 4012 #endif 4013 #if defined(PETSC_HAVE_UMFPACK) 4014 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_umfpack_C",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 4015 #endif 4016 #if defined(PETSC_HAVE_CHOLMOD) 4017 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_cholmod_C",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 4018 #endif 4019 #if defined(PETSC_HAVE_LUSOL) 4020 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_lusol_C",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 4021 #endif 4022 #if defined(PETSC_HAVE_CLIQUE) 4023 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_clique_C",MatGetFactor_aij_clique);CHKERRQ(ierr); 4024 #endif 4025 4026 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 4027 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactorAvailable_petsc_C",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 4028 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_bas_C",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 4029 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4030 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4031 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4032 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4033 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4034 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4035 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4036 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4037 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4038 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4039 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4040 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4041 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4042 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4043 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4044 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4045 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4046 PetscFunctionReturn(0); 4047 } 4048 4049 #undef __FUNCT__ 4050 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 4051 /* 4052 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4053 */ 4054 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4055 { 4056 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4057 PetscErrorCode ierr; 4058 PetscInt i,m = A->rmap->n; 4059 4060 PetscFunctionBegin; 4061 c = (Mat_SeqAIJ*)C->data; 4062 4063 C->factortype = A->factortype; 4064 c->row = 0; 4065 c->col = 0; 4066 c->icol = 0; 4067 c->reallocs = 0; 4068 4069 C->assembled = PETSC_TRUE; 4070 4071 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4072 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4073 4074 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 4075 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4076 for (i=0; i<m; i++) { 4077 c->imax[i] = a->imax[i]; 4078 c->ilen[i] = a->ilen[i]; 4079 } 4080 4081 /* allocate the matrix space */ 4082 if (mallocmatspace) { 4083 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 4084 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4085 4086 c->singlemalloc = PETSC_TRUE; 4087 4088 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4089 if (m > 0) { 4090 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4091 if (cpvalues == MAT_COPY_VALUES) { 4092 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4093 } else { 4094 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4095 } 4096 } 4097 } 4098 4099 c->ignorezeroentries = a->ignorezeroentries; 4100 c->roworiented = a->roworiented; 4101 c->nonew = a->nonew; 4102 if (a->diag) { 4103 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 4104 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4105 for (i=0; i<m; i++) { 4106 c->diag[i] = a->diag[i]; 4107 } 4108 } else c->diag = 0; 4109 4110 c->solve_work = 0; 4111 c->saved_values = 0; 4112 c->idiag = 0; 4113 c->ssor_work = 0; 4114 c->keepnonzeropattern = a->keepnonzeropattern; 4115 c->free_a = PETSC_TRUE; 4116 c->free_ij = PETSC_TRUE; 4117 c->xtoy = 0; 4118 c->XtoY = 0; 4119 4120 c->rmax = a->rmax; 4121 c->nz = a->nz; 4122 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4123 C->preallocated = PETSC_TRUE; 4124 4125 c->compressedrow.use = a->compressedrow.use; 4126 c->compressedrow.nrows = a->compressedrow.nrows; 4127 c->compressedrow.check = a->compressedrow.check; 4128 if (a->compressedrow.use) { 4129 i = a->compressedrow.nrows; 4130 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 4131 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4132 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4133 } else { 4134 c->compressedrow.use = PETSC_FALSE; 4135 c->compressedrow.i = NULL; 4136 c->compressedrow.rindex = NULL; 4137 } 4138 C->same_nonzero = A->same_nonzero; 4139 4140 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4141 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4142 PetscFunctionReturn(0); 4143 } 4144 4145 #undef __FUNCT__ 4146 #define __FUNCT__ "MatDuplicate_SeqAIJ" 4147 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4148 { 4149 PetscErrorCode ierr; 4150 4151 PetscFunctionBegin; 4152 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4153 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4154 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 4155 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4156 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4157 PetscFunctionReturn(0); 4158 } 4159 4160 #undef __FUNCT__ 4161 #define __FUNCT__ "MatLoad_SeqAIJ" 4162 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4163 { 4164 Mat_SeqAIJ *a; 4165 PetscErrorCode ierr; 4166 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4167 int fd; 4168 PetscMPIInt size; 4169 MPI_Comm comm; 4170 PetscInt bs = 1; 4171 4172 PetscFunctionBegin; 4173 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4174 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4175 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4176 4177 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4178 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4179 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4180 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 4181 4182 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4183 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4184 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4185 M = header[1]; N = header[2]; nz = header[3]; 4186 4187 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4188 4189 /* read in row lengths */ 4190 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 4191 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4192 4193 /* check if sum of rowlengths is same as nz */ 4194 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4195 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); 4196 4197 /* set global size if not set already*/ 4198 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4199 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4200 } else { 4201 /* if sizes and type are already set, check if the vector global sizes are correct */ 4202 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4203 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4204 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4205 } 4206 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); 4207 } 4208 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4209 a = (Mat_SeqAIJ*)newMat->data; 4210 4211 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4212 4213 /* read in nonzero values */ 4214 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4215 4216 /* set matrix "i" values */ 4217 a->i[0] = 0; 4218 for (i=1; i<= M; i++) { 4219 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4220 a->ilen[i-1] = rowlengths[i-1]; 4221 } 4222 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4223 4224 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4225 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4226 PetscFunctionReturn(0); 4227 } 4228 4229 #undef __FUNCT__ 4230 #define __FUNCT__ "MatEqual_SeqAIJ" 4231 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4232 { 4233 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4234 PetscErrorCode ierr; 4235 #if defined(PETSC_USE_COMPLEX) 4236 PetscInt k; 4237 #endif 4238 4239 PetscFunctionBegin; 4240 /* If the matrix dimensions are not equal,or no of nonzeros */ 4241 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4242 *flg = PETSC_FALSE; 4243 PetscFunctionReturn(0); 4244 } 4245 4246 /* if the a->i are the same */ 4247 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4248 if (!*flg) PetscFunctionReturn(0); 4249 4250 /* if a->j are the same */ 4251 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4252 if (!*flg) PetscFunctionReturn(0); 4253 4254 /* if a->a are the same */ 4255 #if defined(PETSC_USE_COMPLEX) 4256 for (k=0; k<a->nz; k++) { 4257 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4258 *flg = PETSC_FALSE; 4259 PetscFunctionReturn(0); 4260 } 4261 } 4262 #else 4263 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4264 #endif 4265 PetscFunctionReturn(0); 4266 } 4267 4268 #undef __FUNCT__ 4269 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4270 /*@ 4271 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4272 provided by the user. 4273 4274 Collective on MPI_Comm 4275 4276 Input Parameters: 4277 + comm - must be an MPI communicator of size 1 4278 . m - number of rows 4279 . n - number of columns 4280 . i - row indices 4281 . j - column indices 4282 - a - matrix values 4283 4284 Output Parameter: 4285 . mat - the matrix 4286 4287 Level: intermediate 4288 4289 Notes: 4290 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4291 once the matrix is destroyed and not before 4292 4293 You cannot set new nonzero locations into this matrix, that will generate an error. 4294 4295 The i and j indices are 0 based 4296 4297 The format which is used for the sparse matrix input, is equivalent to a 4298 row-major ordering.. i.e for the following matrix, the input data expected is 4299 as shown: 4300 4301 1 0 0 4302 2 0 3 4303 4 5 6 4304 4305 i = {0,1,3,6} [size = nrow+1 = 3+1] 4306 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4307 v = {1,2,3,4,5,6} [size = nz = 6] 4308 4309 4310 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4311 4312 @*/ 4313 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 4314 { 4315 PetscErrorCode ierr; 4316 PetscInt ii; 4317 Mat_SeqAIJ *aij; 4318 #if defined(PETSC_USE_DEBUG) 4319 PetscInt jj; 4320 #endif 4321 4322 PetscFunctionBegin; 4323 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4324 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4325 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4326 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4327 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4328 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4329 aij = (Mat_SeqAIJ*)(*mat)->data; 4330 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 4331 4332 aij->i = i; 4333 aij->j = j; 4334 aij->a = a; 4335 aij->singlemalloc = PETSC_FALSE; 4336 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4337 aij->free_a = PETSC_FALSE; 4338 aij->free_ij = PETSC_FALSE; 4339 4340 for (ii=0; ii<m; ii++) { 4341 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4342 #if defined(PETSC_USE_DEBUG) 4343 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]); 4344 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4345 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); 4346 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); 4347 } 4348 #endif 4349 } 4350 #if defined(PETSC_USE_DEBUG) 4351 for (ii=0; ii<aij->i[m]; ii++) { 4352 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 4353 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]); 4354 } 4355 #endif 4356 4357 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4358 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4359 PetscFunctionReturn(0); 4360 } 4361 #undef __FUNCT__ 4362 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4363 /*@C 4364 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4365 provided by the user. 4366 4367 Collective on MPI_Comm 4368 4369 Input Parameters: 4370 + comm - must be an MPI communicator of size 1 4371 . m - number of rows 4372 . n - number of columns 4373 . i - row indices 4374 . j - column indices 4375 . a - matrix values 4376 . nz - number of nonzeros 4377 - idx - 0 or 1 based 4378 4379 Output Parameter: 4380 . mat - the matrix 4381 4382 Level: intermediate 4383 4384 Notes: 4385 The i and j indices are 0 based 4386 4387 The format which is used for the sparse matrix input, is equivalent to a 4388 row-major ordering.. i.e for the following matrix, the input data expected is 4389 as shown: 4390 4391 1 0 0 4392 2 0 3 4393 4 5 6 4394 4395 i = {0,1,1,2,2,2} 4396 j = {0,0,2,0,1,2} 4397 v = {1,2,3,4,5,6} 4398 4399 4400 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4401 4402 @*/ 4403 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4404 { 4405 PetscErrorCode ierr; 4406 PetscInt ii, *nnz, one = 1,row,col; 4407 4408 4409 PetscFunctionBegin; 4410 ierr = PetscMalloc(m*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 4411 ierr = PetscMemzero(nnz,m*sizeof(PetscInt));CHKERRQ(ierr); 4412 for (ii = 0; ii < nz; ii++) { 4413 nnz[i[ii]] += 1; 4414 } 4415 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4416 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4417 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4418 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4419 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4420 for (ii = 0; ii < nz; ii++) { 4421 if (idx) { 4422 row = i[ii] - 1; 4423 col = j[ii] - 1; 4424 } else { 4425 row = i[ii]; 4426 col = j[ii]; 4427 } 4428 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4429 } 4430 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4431 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4432 ierr = PetscFree(nnz);CHKERRQ(ierr); 4433 PetscFunctionReturn(0); 4434 } 4435 4436 #undef __FUNCT__ 4437 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4438 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4439 { 4440 PetscErrorCode ierr; 4441 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4442 4443 PetscFunctionBegin; 4444 if (coloring->ctype == IS_COLORING_GLOBAL) { 4445 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4446 a->coloring = coloring; 4447 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4448 PetscInt i,*larray; 4449 ISColoring ocoloring; 4450 ISColoringValue *colors; 4451 4452 /* set coloring for diagonal portion */ 4453 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 4454 for (i=0; i<A->cmap->n; i++) larray[i] = i; 4455 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,NULL,larray);CHKERRQ(ierr); 4456 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 4457 for (i=0; i<A->cmap->n; i++) colors[i] = coloring->colors[larray[i]]; 4458 ierr = PetscFree(larray);CHKERRQ(ierr); 4459 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4460 a->coloring = ocoloring; 4461 } 4462 PetscFunctionReturn(0); 4463 } 4464 4465 #undef __FUNCT__ 4466 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4467 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4468 { 4469 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4470 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4471 MatScalar *v = a->a; 4472 PetscScalar *values = (PetscScalar*)advalues; 4473 ISColoringValue *color; 4474 4475 PetscFunctionBegin; 4476 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4477 color = a->coloring->colors; 4478 /* loop over rows */ 4479 for (i=0; i<m; i++) { 4480 nz = ii[i+1] - ii[i]; 4481 /* loop over columns putting computed value into matrix */ 4482 for (j=0; j<nz; j++) *v++ = values[color[*jj++]]; 4483 values += nl; /* jump to next row of derivatives */ 4484 } 4485 PetscFunctionReturn(0); 4486 } 4487 4488 #undef __FUNCT__ 4489 #define __FUNCT__ "MatSeqAIJInvalidateDiagonal" 4490 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4491 { 4492 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4493 PetscErrorCode ierr; 4494 4495 PetscFunctionBegin; 4496 a->idiagvalid = PETSC_FALSE; 4497 a->ibdiagvalid = PETSC_FALSE; 4498 4499 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4500 PetscFunctionReturn(0); 4501 } 4502 4503 /* 4504 Special version for direct calls from Fortran 4505 */ 4506 #include <petsc-private/fortranimpl.h> 4507 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4508 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4509 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4510 #define matsetvaluesseqaij_ matsetvaluesseqaij 4511 #endif 4512 4513 /* Change these macros so can be used in void function */ 4514 #undef CHKERRQ 4515 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4516 #undef SETERRQ2 4517 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4518 #undef SETERRQ3 4519 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4520 4521 #undef __FUNCT__ 4522 #define __FUNCT__ "matsetvaluesseqaij_" 4523 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) 4524 { 4525 Mat A = *AA; 4526 PetscInt m = *mm, n = *nn; 4527 InsertMode is = *isis; 4528 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4529 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4530 PetscInt *imax,*ai,*ailen; 4531 PetscErrorCode ierr; 4532 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4533 MatScalar *ap,value,*aa; 4534 PetscBool ignorezeroentries = a->ignorezeroentries; 4535 PetscBool roworiented = a->roworiented; 4536 4537 PetscFunctionBegin; 4538 MatCheckPreallocated(A,1); 4539 imax = a->imax; 4540 ai = a->i; 4541 ailen = a->ilen; 4542 aj = a->j; 4543 aa = a->a; 4544 4545 for (k=0; k<m; k++) { /* loop over added rows */ 4546 row = im[k]; 4547 if (row < 0) continue; 4548 #if defined(PETSC_USE_DEBUG) 4549 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4550 #endif 4551 rp = aj + ai[row]; ap = aa + ai[row]; 4552 rmax = imax[row]; nrow = ailen[row]; 4553 low = 0; 4554 high = nrow; 4555 for (l=0; l<n; l++) { /* loop over added columns */ 4556 if (in[l] < 0) continue; 4557 #if defined(PETSC_USE_DEBUG) 4558 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4559 #endif 4560 col = in[l]; 4561 if (roworiented) value = v[l + k*n]; 4562 else value = v[k + l*m]; 4563 4564 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4565 4566 if (col <= lastcol) low = 0; 4567 else high = nrow; 4568 lastcol = col; 4569 while (high-low > 5) { 4570 t = (low+high)/2; 4571 if (rp[t] > col) high = t; 4572 else low = t; 4573 } 4574 for (i=low; i<high; i++) { 4575 if (rp[i] > col) break; 4576 if (rp[i] == col) { 4577 if (is == ADD_VALUES) ap[i] += value; 4578 else ap[i] = value; 4579 goto noinsert; 4580 } 4581 } 4582 if (value == 0.0 && ignorezeroentries) goto noinsert; 4583 if (nonew == 1) goto noinsert; 4584 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4585 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4586 N = nrow++ - 1; a->nz++; high++; 4587 /* shift up all the later entries in this row */ 4588 for (ii=N; ii>=i; ii--) { 4589 rp[ii+1] = rp[ii]; 4590 ap[ii+1] = ap[ii]; 4591 } 4592 rp[i] = col; 4593 ap[i] = value; 4594 noinsert:; 4595 low = i + 1; 4596 } 4597 ailen[row] = nrow; 4598 } 4599 A->same_nonzero = PETSC_FALSE; 4600 PetscFunctionReturnVoid(); 4601 } 4602 4603 4604