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