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