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