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