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