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