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