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 for (i=0; i<m; i++) { 1498 len = ai[i+1]-ai[i]; 1499 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1500 array += len; 1501 aj += len; 1502 } 1503 } else { 1504 C = *B; 1505 } 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 PetscInt *irow,*icol,nrows,ncols; 1699 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1700 MatScalar *a_new,*mat_a; 1701 Mat C; 1702 PetscTruth stride; 1703 1704 PetscFunctionBegin; 1705 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1706 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1707 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1708 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1709 1710 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1711 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1712 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1713 1714 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1715 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1716 if (stride && step == 1) { 1717 /* special case of contiguous rows */ 1718 ierr = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1719 starts = lens + nrows; 1720 /* loop over new rows determining lens and starting points */ 1721 for (i=0; i<nrows; i++) { 1722 kstart = ai[irow[i]]; 1723 kend = kstart + ailen[irow[i]]; 1724 for (k=kstart; k<kend; k++) { 1725 if (aj[k] >= first) { 1726 starts[i] = k; 1727 break; 1728 } 1729 } 1730 sum = 0; 1731 while (k < kend) { 1732 if (aj[k++] >= first+ncols) break; 1733 sum++; 1734 } 1735 lens[i] = sum; 1736 } 1737 /* create submatrix */ 1738 if (scall == MAT_REUSE_MATRIX) { 1739 PetscInt n_cols,n_rows; 1740 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1741 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1742 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1743 C = *B; 1744 } else { 1745 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1746 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1747 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1748 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1749 } 1750 c = (Mat_SeqAIJ*)C->data; 1751 1752 /* loop over rows inserting into submatrix */ 1753 a_new = c->a; 1754 j_new = c->j; 1755 i_new = c->i; 1756 1757 for (i=0; i<nrows; i++) { 1758 ii = starts[i]; 1759 lensi = lens[i]; 1760 for (k=0; k<lensi; k++) { 1761 *j_new++ = aj[ii+k] - first; 1762 } 1763 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1764 a_new += lensi; 1765 i_new[i+1] = i_new[i] + lensi; 1766 c->ilen[i] = lensi; 1767 } 1768 ierr = PetscFree(lens);CHKERRQ(ierr); 1769 } else { 1770 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1771 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1772 1773 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1774 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1775 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1776 /* determine lens of each row */ 1777 for (i=0; i<nrows; i++) { 1778 kstart = ai[irow[i]]; 1779 kend = kstart + a->ilen[irow[i]]; 1780 lens[i] = 0; 1781 for (k=kstart; k<kend; k++) { 1782 if (smap[aj[k]]) { 1783 lens[i]++; 1784 } 1785 } 1786 } 1787 /* Create and fill new matrix */ 1788 if (scall == MAT_REUSE_MATRIX) { 1789 PetscTruth equal; 1790 1791 c = (Mat_SeqAIJ *)((*B)->data); 1792 if ((*B)->rmap.n != nrows || (*B)->cmap.n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1793 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap.n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1794 if (!equal) { 1795 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1796 } 1797 ierr = PetscMemzero(c->ilen,(*B)->rmap.n*sizeof(PetscInt));CHKERRQ(ierr); 1798 C = *B; 1799 } else { 1800 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1801 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1802 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1803 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1804 } 1805 c = (Mat_SeqAIJ *)(C->data); 1806 for (i=0; i<nrows; i++) { 1807 row = irow[i]; 1808 kstart = ai[row]; 1809 kend = kstart + a->ilen[row]; 1810 mat_i = c->i[i]; 1811 mat_j = c->j + mat_i; 1812 mat_a = c->a + mat_i; 1813 mat_ilen = c->ilen + i; 1814 for (k=kstart; k<kend; k++) { 1815 if ((tcol=smap[a->j[k]])) { 1816 *mat_j++ = tcol - 1; 1817 *mat_a++ = a->a[k]; 1818 (*mat_ilen)++; 1819 1820 } 1821 } 1822 } 1823 /* Free work space */ 1824 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1825 ierr = PetscFree(smap);CHKERRQ(ierr); 1826 ierr = PetscFree(lens);CHKERRQ(ierr); 1827 } 1828 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1829 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1830 1831 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1832 *B = C; 1833 PetscFunctionReturn(0); 1834 } 1835 1836 /* 1837 */ 1838 #undef __FUNCT__ 1839 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1840 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1841 { 1842 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1843 PetscErrorCode ierr; 1844 Mat outA; 1845 PetscTruth row_identity,col_identity; 1846 1847 PetscFunctionBegin; 1848 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1849 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1850 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1851 1852 outA = inA; 1853 inA->factor = MAT_FACTOR_LU; 1854 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1855 if (a->row) { ierr = ISDestroy(a->row); CHKERRQ(ierr);} 1856 a->row = row; 1857 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1858 if (a->col) { ierr = ISDestroy(a->col); CHKERRQ(ierr);} 1859 a->col = col; 1860 1861 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1862 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1863 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1864 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1865 1866 if (!a->solve_work) { /* this matrix may have been factored before */ 1867 ierr = PetscMalloc((inA->rmap.n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1868 ierr = PetscLogObjectMemory(inA, (inA->rmap.n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1869 } 1870 1871 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1872 if (row_identity && col_identity) { 1873 ierr = MatLUFactorNumeric_SeqAIJ(inA,info,&outA);CHKERRQ(ierr); 1874 } else { 1875 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(inA,info,&outA);CHKERRQ(ierr); 1876 } 1877 PetscFunctionReturn(0); 1878 } 1879 1880 #include "petscblaslapack.h" 1881 #undef __FUNCT__ 1882 #define __FUNCT__ "MatScale_SeqAIJ" 1883 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1884 { 1885 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1886 PetscScalar oalpha = alpha; 1887 PetscErrorCode ierr; 1888 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1889 1890 PetscFunctionBegin; 1891 BLASscal_(&bnz,&oalpha,a->a,&one); 1892 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1893 PetscFunctionReturn(0); 1894 } 1895 1896 #undef __FUNCT__ 1897 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1898 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1899 { 1900 PetscErrorCode ierr; 1901 PetscInt i; 1902 1903 PetscFunctionBegin; 1904 if (scall == MAT_INITIAL_MATRIX) { 1905 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1906 } 1907 1908 for (i=0; i<n; i++) { 1909 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1910 } 1911 PetscFunctionReturn(0); 1912 } 1913 1914 #undef __FUNCT__ 1915 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1916 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1917 { 1918 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1919 PetscErrorCode ierr; 1920 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val; 1921 PetscInt start,end,*ai,*aj; 1922 PetscBT table; 1923 1924 PetscFunctionBegin; 1925 m = A->rmap.n; 1926 ai = a->i; 1927 aj = a->j; 1928 1929 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1930 1931 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1932 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1933 1934 for (i=0; i<is_max; i++) { 1935 /* Initialize the two local arrays */ 1936 isz = 0; 1937 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1938 1939 /* Extract the indices, assume there can be duplicate entries */ 1940 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1941 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1942 1943 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1944 for (j=0; j<n ; ++j){ 1945 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1946 } 1947 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1948 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1949 1950 k = 0; 1951 for (j=0; j<ov; j++){ /* for each overlap */ 1952 n = isz; 1953 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1954 row = nidx[k]; 1955 start = ai[row]; 1956 end = ai[row+1]; 1957 for (l = start; l<end ; l++){ 1958 val = aj[l] ; 1959 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1960 } 1961 } 1962 } 1963 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1964 } 1965 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1966 ierr = PetscFree(nidx);CHKERRQ(ierr); 1967 PetscFunctionReturn(0); 1968 } 1969 1970 /* -------------------------------------------------------------- */ 1971 #undef __FUNCT__ 1972 #define __FUNCT__ "MatPermute_SeqAIJ" 1973 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1974 { 1975 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1976 PetscErrorCode ierr; 1977 PetscInt i,nz,m = A->rmap.n,n = A->cmap.n,*col; 1978 PetscInt *row,*cnew,j,*lens; 1979 IS icolp,irowp; 1980 PetscInt *cwork; 1981 PetscScalar *vwork; 1982 1983 PetscFunctionBegin; 1984 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1985 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1986 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1987 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1988 1989 /* determine lengths of permuted rows */ 1990 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1991 for (i=0; i<m; i++) { 1992 lens[row[i]] = a->i[i+1] - a->i[i]; 1993 } 1994 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1995 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 1996 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 1997 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 1998 ierr = PetscFree(lens);CHKERRQ(ierr); 1999 2000 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2001 for (i=0; i<m; i++) { 2002 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2003 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2004 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2005 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2006 } 2007 ierr = PetscFree(cnew);CHKERRQ(ierr); 2008 (*B)->assembled = PETSC_FALSE; 2009 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2010 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2011 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2012 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2013 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2014 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2015 PetscFunctionReturn(0); 2016 } 2017 2018 #undef __FUNCT__ 2019 #define __FUNCT__ "MatCopy_SeqAIJ" 2020 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2021 { 2022 PetscErrorCode ierr; 2023 2024 PetscFunctionBegin; 2025 /* If the two matrices have the same copy implementation, use fast copy. */ 2026 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2027 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2028 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2029 2030 if (a->i[A->rmap.n] != b->i[B->rmap.n]) { 2031 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2032 } 2033 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap.n])*sizeof(PetscScalar));CHKERRQ(ierr); 2034 } else { 2035 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2036 } 2037 PetscFunctionReturn(0); 2038 } 2039 2040 #undef __FUNCT__ 2041 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2042 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2043 { 2044 PetscErrorCode ierr; 2045 2046 PetscFunctionBegin; 2047 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2048 PetscFunctionReturn(0); 2049 } 2050 2051 #undef __FUNCT__ 2052 #define __FUNCT__ "MatGetArray_SeqAIJ" 2053 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2054 { 2055 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2056 PetscFunctionBegin; 2057 *array = a->a; 2058 PetscFunctionReturn(0); 2059 } 2060 2061 #undef __FUNCT__ 2062 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2063 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2064 { 2065 PetscFunctionBegin; 2066 PetscFunctionReturn(0); 2067 } 2068 2069 #undef __FUNCT__ 2070 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2071 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2072 { 2073 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2074 PetscErrorCode ierr; 2075 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2076 PetscScalar dx,*y,*xx,*w3_array; 2077 PetscScalar *vscale_array; 2078 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2079 Vec w1,w2,w3; 2080 void *fctx = coloring->fctx; 2081 PetscTruth flg; 2082 2083 PetscFunctionBegin; 2084 if (!coloring->w1) { 2085 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2086 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2087 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2088 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2089 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2090 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2091 } 2092 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2093 2094 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2095 ierr = PetscOptionsHasName(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 2096 if (flg) { 2097 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2098 } else { 2099 PetscTruth assembled; 2100 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2101 if (assembled) { 2102 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2103 } 2104 } 2105 2106 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2107 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2108 2109 /* 2110 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2111 coloring->F for the coarser grids from the finest 2112 */ 2113 if (coloring->F) { 2114 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2115 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2116 if (m1 != m2) { 2117 coloring->F = 0; 2118 } 2119 } 2120 2121 if (coloring->F) { 2122 w1 = coloring->F; 2123 coloring->F = 0; 2124 } else { 2125 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2126 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2127 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2128 } 2129 2130 /* 2131 Compute all the scale factors and share with other processors 2132 */ 2133 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2134 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2135 for (k=0; k<coloring->ncolors; k++) { 2136 /* 2137 Loop over each column associated with color adding the 2138 perturbation to the vector w3. 2139 */ 2140 for (l=0; l<coloring->ncolumns[k]; l++) { 2141 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2142 dx = xx[col]; 2143 if (dx == 0.0) dx = 1.0; 2144 #if !defined(PETSC_USE_COMPLEX) 2145 if (dx < umin && dx >= 0.0) dx = umin; 2146 else if (dx < 0.0 && dx > -umin) dx = -umin; 2147 #else 2148 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2149 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2150 #endif 2151 dx *= epsilon; 2152 vscale_array[col] = 1.0/dx; 2153 } 2154 } 2155 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2156 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2157 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2158 2159 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2160 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2161 2162 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2163 else vscaleforrow = coloring->columnsforrow; 2164 2165 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2166 /* 2167 Loop over each color 2168 */ 2169 for (k=0; k<coloring->ncolors; k++) { 2170 coloring->currentcolor = k; 2171 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2172 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2173 /* 2174 Loop over each column associated with color adding the 2175 perturbation to the vector w3. 2176 */ 2177 for (l=0; l<coloring->ncolumns[k]; l++) { 2178 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2179 dx = xx[col]; 2180 if (dx == 0.0) dx = 1.0; 2181 #if !defined(PETSC_USE_COMPLEX) 2182 if (dx < umin && dx >= 0.0) dx = umin; 2183 else if (dx < 0.0 && dx > -umin) dx = -umin; 2184 #else 2185 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2186 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2187 #endif 2188 dx *= epsilon; 2189 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2190 w3_array[col] += dx; 2191 } 2192 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2193 2194 /* 2195 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2196 */ 2197 2198 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2199 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2200 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2201 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2202 2203 /* 2204 Loop over rows of vector, putting results into Jacobian matrix 2205 */ 2206 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2207 for (l=0; l<coloring->nrows[k]; l++) { 2208 row = coloring->rows[k][l]; 2209 col = coloring->columnsforrow[k][l]; 2210 y[row] *= vscale_array[vscaleforrow[k][l]]; 2211 srow = row + start; 2212 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2213 } 2214 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2215 } 2216 coloring->currentcolor = k; 2217 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2218 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2219 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2220 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2221 PetscFunctionReturn(0); 2222 } 2223 2224 #include "petscblaslapack.h" 2225 #undef __FUNCT__ 2226 #define __FUNCT__ "MatAXPY_SeqAIJ" 2227 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2228 { 2229 PetscErrorCode ierr; 2230 PetscInt i; 2231 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2232 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2233 2234 PetscFunctionBegin; 2235 if (str == SAME_NONZERO_PATTERN) { 2236 PetscScalar alpha = a; 2237 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2238 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2239 if (y->xtoy && y->XtoY != X) { 2240 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2241 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2242 } 2243 if (!y->xtoy) { /* get xtoy */ 2244 ierr = MatAXPYGetxtoy_Private(X->rmap.n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2245 y->XtoY = X; 2246 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2247 } 2248 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2249 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2250 } else { 2251 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2252 } 2253 PetscFunctionReturn(0); 2254 } 2255 2256 #undef __FUNCT__ 2257 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2258 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2259 { 2260 PetscFunctionBegin; 2261 PetscFunctionReturn(0); 2262 } 2263 2264 #undef __FUNCT__ 2265 #define __FUNCT__ "MatConjugate_SeqAIJ" 2266 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2267 { 2268 #if defined(PETSC_USE_COMPLEX) 2269 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2270 PetscInt i,nz; 2271 PetscScalar *a; 2272 2273 PetscFunctionBegin; 2274 nz = aij->nz; 2275 a = aij->a; 2276 for (i=0; i<nz; i++) { 2277 a[i] = PetscConj(a[i]); 2278 } 2279 #else 2280 PetscFunctionBegin; 2281 #endif 2282 PetscFunctionReturn(0); 2283 } 2284 2285 #undef __FUNCT__ 2286 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2287 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2288 { 2289 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2290 PetscErrorCode ierr; 2291 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2292 PetscReal atmp; 2293 PetscScalar *x; 2294 MatScalar *aa; 2295 2296 PetscFunctionBegin; 2297 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2298 aa = a->a; 2299 ai = a->i; 2300 aj = a->j; 2301 2302 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2303 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2304 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2305 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2306 for (i=0; i<m; i++) { 2307 ncols = ai[1] - ai[0]; ai++; 2308 x[i] = 0.0; 2309 for (j=0; j<ncols; j++){ 2310 atmp = PetscAbsScalar(*aa); 2311 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2312 aa++; aj++; 2313 } 2314 } 2315 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2316 PetscFunctionReturn(0); 2317 } 2318 2319 #undef __FUNCT__ 2320 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2321 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2322 { 2323 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2324 PetscErrorCode ierr; 2325 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2326 PetscScalar *x; 2327 MatScalar *aa; 2328 2329 PetscFunctionBegin; 2330 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2331 aa = a->a; 2332 ai = a->i; 2333 aj = a->j; 2334 2335 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2336 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2337 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2338 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2339 for (i=0; i<m; i++) { 2340 ncols = ai[1] - ai[0]; ai++; 2341 if (ncols == A->cmap.n) { /* row is dense */ 2342 x[i] = *aa; if (idx) idx[i] = 0; 2343 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2344 x[i] = 0.0; 2345 if (idx) { 2346 idx[i] = 0; /* in case ncols is zero */ 2347 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2348 if (aj[j] > j) { 2349 idx[i] = j; 2350 break; 2351 } 2352 } 2353 } 2354 } 2355 for (j=0; j<ncols; j++){ 2356 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2357 aa++; aj++; 2358 } 2359 } 2360 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2361 PetscFunctionReturn(0); 2362 } 2363 2364 #undef __FUNCT__ 2365 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2366 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2367 { 2368 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2369 PetscErrorCode ierr; 2370 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2371 PetscScalar *x; 2372 MatScalar *aa; 2373 2374 PetscFunctionBegin; 2375 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2376 aa = a->a; 2377 ai = a->i; 2378 aj = a->j; 2379 2380 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2381 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2382 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2383 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2384 for (i=0; i<m; i++) { 2385 ncols = ai[1] - ai[0]; ai++; 2386 if (ncols == A->cmap.n) { /* row is dense */ 2387 x[i] = *aa; if (idx) idx[i] = 0; 2388 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2389 x[i] = 0.0; 2390 if (idx) { /* find first implicit 0.0 in the row */ 2391 idx[i] = 0; /* in case ncols is zero */ 2392 for (j=0;j<ncols;j++) { 2393 if (aj[j] > j) { 2394 idx[i] = j; 2395 break; 2396 } 2397 } 2398 } 2399 } 2400 for (j=0; j<ncols; j++){ 2401 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2402 aa++; aj++; 2403 } 2404 } 2405 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2406 PetscFunctionReturn(0); 2407 } 2408 2409 /* -------------------------------------------------------------------*/ 2410 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2411 MatGetRow_SeqAIJ, 2412 MatRestoreRow_SeqAIJ, 2413 MatMult_SeqAIJ, 2414 /* 4*/ MatMultAdd_SeqAIJ, 2415 MatMultTranspose_SeqAIJ, 2416 MatMultTransposeAdd_SeqAIJ, 2417 MatSolve_SeqAIJ, 2418 MatSolveAdd_SeqAIJ, 2419 MatSolveTranspose_SeqAIJ, 2420 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2421 MatLUFactor_SeqAIJ, 2422 0, 2423 MatRelax_SeqAIJ, 2424 MatTranspose_SeqAIJ, 2425 /*15*/ MatGetInfo_SeqAIJ, 2426 MatEqual_SeqAIJ, 2427 MatGetDiagonal_SeqAIJ, 2428 MatDiagonalScale_SeqAIJ, 2429 MatNorm_SeqAIJ, 2430 /*20*/ 0, 2431 MatAssemblyEnd_SeqAIJ, 2432 MatCompress_SeqAIJ, 2433 MatSetOption_SeqAIJ, 2434 MatZeroEntries_SeqAIJ, 2435 /*25*/ MatZeroRows_SeqAIJ, 2436 MatLUFactorSymbolic_SeqAIJ, 2437 MatLUFactorNumeric_SeqAIJ, 2438 MatCholeskyFactorSymbolic_SeqAIJ, 2439 MatCholeskyFactorNumeric_SeqAIJ, 2440 /*30*/ MatSetUpPreallocation_SeqAIJ, 2441 MatILUFactorSymbolic_SeqAIJ, 2442 MatICCFactorSymbolic_SeqAIJ, 2443 MatGetArray_SeqAIJ, 2444 MatRestoreArray_SeqAIJ, 2445 /*35*/ MatDuplicate_SeqAIJ, 2446 0, 2447 0, 2448 MatILUFactor_SeqAIJ, 2449 0, 2450 /*40*/ MatAXPY_SeqAIJ, 2451 MatGetSubMatrices_SeqAIJ, 2452 MatIncreaseOverlap_SeqAIJ, 2453 MatGetValues_SeqAIJ, 2454 MatCopy_SeqAIJ, 2455 /*45*/ MatGetRowMax_SeqAIJ, 2456 MatScale_SeqAIJ, 2457 0, 2458 MatDiagonalSet_SeqAIJ, 2459 MatILUDTFactor_SeqAIJ, 2460 /*50*/ MatSetBlockSize_SeqAIJ, 2461 MatGetRowIJ_SeqAIJ, 2462 MatRestoreRowIJ_SeqAIJ, 2463 MatGetColumnIJ_SeqAIJ, 2464 MatRestoreColumnIJ_SeqAIJ, 2465 /*55*/ MatFDColoringCreate_SeqAIJ, 2466 0, 2467 0, 2468 MatPermute_SeqAIJ, 2469 0, 2470 /*60*/ 0, 2471 MatDestroy_SeqAIJ, 2472 MatView_SeqAIJ, 2473 0, 2474 0, 2475 /*65*/ 0, 2476 0, 2477 0, 2478 0, 2479 0, 2480 /*70*/ MatGetRowMaxAbs_SeqAIJ, 2481 0, 2482 MatSetColoring_SeqAIJ, 2483 #if defined(PETSC_HAVE_ADIC) 2484 MatSetValuesAdic_SeqAIJ, 2485 #else 2486 0, 2487 #endif 2488 MatSetValuesAdifor_SeqAIJ, 2489 /*75*/ 0, 2490 0, 2491 0, 2492 0, 2493 0, 2494 /*80*/ 0, 2495 0, 2496 0, 2497 0, 2498 MatLoad_SeqAIJ, 2499 /*85*/ MatIsSymmetric_SeqAIJ, 2500 MatIsHermitian_SeqAIJ, 2501 0, 2502 0, 2503 0, 2504 /*90*/ MatMatMult_SeqAIJ_SeqAIJ, 2505 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2506 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2507 MatPtAP_Basic, 2508 MatPtAPSymbolic_SeqAIJ, 2509 /*95*/ MatPtAPNumeric_SeqAIJ, 2510 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2511 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2512 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2513 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2514 /*100*/MatPtAPNumeric_SeqAIJ_SeqAIJ, 2515 0, 2516 0, 2517 MatConjugate_SeqAIJ, 2518 0, 2519 /*105*/MatSetValuesRow_SeqAIJ, 2520 MatRealPart_SeqAIJ, 2521 MatImaginaryPart_SeqAIJ, 2522 0, 2523 0, 2524 /*110*/MatMatSolve_SeqAIJ, 2525 0, 2526 MatGetRowMin_SeqAIJ, 2527 0, 2528 MatMissingDiagonal_SeqAIJ 2529 }; 2530 2531 EXTERN_C_BEGIN 2532 #undef __FUNCT__ 2533 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2534 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2535 { 2536 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2537 PetscInt i,nz,n; 2538 2539 PetscFunctionBegin; 2540 2541 nz = aij->maxnz; 2542 n = mat->rmap.n; 2543 for (i=0; i<nz; i++) { 2544 aij->j[i] = indices[i]; 2545 } 2546 aij->nz = nz; 2547 for (i=0; i<n; i++) { 2548 aij->ilen[i] = aij->imax[i]; 2549 } 2550 2551 PetscFunctionReturn(0); 2552 } 2553 EXTERN_C_END 2554 2555 #undef __FUNCT__ 2556 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2557 /*@ 2558 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2559 in the matrix. 2560 2561 Input Parameters: 2562 + mat - the SeqAIJ matrix 2563 - indices - the column indices 2564 2565 Level: advanced 2566 2567 Notes: 2568 This can be called if you have precomputed the nonzero structure of the 2569 matrix and want to provide it to the matrix object to improve the performance 2570 of the MatSetValues() operation. 2571 2572 You MUST have set the correct numbers of nonzeros per row in the call to 2573 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2574 2575 MUST be called before any calls to MatSetValues(); 2576 2577 The indices should start with zero, not one. 2578 2579 @*/ 2580 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2581 { 2582 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2583 2584 PetscFunctionBegin; 2585 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2586 PetscValidPointer(indices,2); 2587 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2588 if (f) { 2589 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2590 } else { 2591 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2592 } 2593 PetscFunctionReturn(0); 2594 } 2595 2596 /* ----------------------------------------------------------------------------------------*/ 2597 2598 EXTERN_C_BEGIN 2599 #undef __FUNCT__ 2600 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2601 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2602 { 2603 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2604 PetscErrorCode ierr; 2605 size_t nz = aij->i[mat->rmap.n]; 2606 2607 PetscFunctionBegin; 2608 if (aij->nonew != 1) { 2609 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2610 } 2611 2612 /* allocate space for values if not already there */ 2613 if (!aij->saved_values) { 2614 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2615 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2616 } 2617 2618 /* copy values over */ 2619 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2620 PetscFunctionReturn(0); 2621 } 2622 EXTERN_C_END 2623 2624 #undef __FUNCT__ 2625 #define __FUNCT__ "MatStoreValues" 2626 /*@ 2627 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2628 example, reuse of the linear part of a Jacobian, while recomputing the 2629 nonlinear portion. 2630 2631 Collect on Mat 2632 2633 Input Parameters: 2634 . mat - the matrix (currently only AIJ matrices support this option) 2635 2636 Level: advanced 2637 2638 Common Usage, with SNESSolve(): 2639 $ Create Jacobian matrix 2640 $ Set linear terms into matrix 2641 $ Apply boundary conditions to matrix, at this time matrix must have 2642 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2643 $ boundary conditions again will not change the nonzero structure 2644 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2645 $ ierr = MatStoreValues(mat); 2646 $ Call SNESSetJacobian() with matrix 2647 $ In your Jacobian routine 2648 $ ierr = MatRetrieveValues(mat); 2649 $ Set nonlinear terms in matrix 2650 2651 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2652 $ // build linear portion of Jacobian 2653 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2654 $ ierr = MatStoreValues(mat); 2655 $ loop over nonlinear iterations 2656 $ ierr = MatRetrieveValues(mat); 2657 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2658 $ // call MatAssemblyBegin/End() on matrix 2659 $ Solve linear system with Jacobian 2660 $ endloop 2661 2662 Notes: 2663 Matrix must already be assemblied before calling this routine 2664 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2665 calling this routine. 2666 2667 When this is called multiple times it overwrites the previous set of stored values 2668 and does not allocated additional space. 2669 2670 .seealso: MatRetrieveValues() 2671 2672 @*/ 2673 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2674 { 2675 PetscErrorCode ierr,(*f)(Mat); 2676 2677 PetscFunctionBegin; 2678 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2679 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2680 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2681 2682 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2683 if (f) { 2684 ierr = (*f)(mat);CHKERRQ(ierr); 2685 } else { 2686 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2687 } 2688 PetscFunctionReturn(0); 2689 } 2690 2691 EXTERN_C_BEGIN 2692 #undef __FUNCT__ 2693 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2694 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2695 { 2696 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2697 PetscErrorCode ierr; 2698 PetscInt nz = aij->i[mat->rmap.n]; 2699 2700 PetscFunctionBegin; 2701 if (aij->nonew != 1) { 2702 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2703 } 2704 if (!aij->saved_values) { 2705 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2706 } 2707 /* copy values over */ 2708 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2709 PetscFunctionReturn(0); 2710 } 2711 EXTERN_C_END 2712 2713 #undef __FUNCT__ 2714 #define __FUNCT__ "MatRetrieveValues" 2715 /*@ 2716 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2717 example, reuse of the linear part of a Jacobian, while recomputing the 2718 nonlinear portion. 2719 2720 Collect on Mat 2721 2722 Input Parameters: 2723 . mat - the matrix (currently on AIJ matrices support this option) 2724 2725 Level: advanced 2726 2727 .seealso: MatStoreValues() 2728 2729 @*/ 2730 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2731 { 2732 PetscErrorCode ierr,(*f)(Mat); 2733 2734 PetscFunctionBegin; 2735 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2736 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2737 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2738 2739 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2740 if (f) { 2741 ierr = (*f)(mat);CHKERRQ(ierr); 2742 } else { 2743 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2744 } 2745 PetscFunctionReturn(0); 2746 } 2747 2748 2749 /* --------------------------------------------------------------------------------*/ 2750 #undef __FUNCT__ 2751 #define __FUNCT__ "MatCreateSeqAIJ" 2752 /*@C 2753 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2754 (the default parallel PETSc format). For good matrix assembly performance 2755 the user should preallocate the matrix storage by setting the parameter nz 2756 (or the array nnz). By setting these parameters accurately, performance 2757 during matrix assembly can be increased by more than a factor of 50. 2758 2759 Collective on MPI_Comm 2760 2761 Input Parameters: 2762 + comm - MPI communicator, set to PETSC_COMM_SELF 2763 . m - number of rows 2764 . n - number of columns 2765 . nz - number of nonzeros per row (same for all rows) 2766 - nnz - array containing the number of nonzeros in the various rows 2767 (possibly different for each row) or PETSC_NULL 2768 2769 Output Parameter: 2770 . A - the matrix 2771 2772 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2773 MatXXXXSetPreallocation() paradgm instead of this routine directly. This is definitely 2774 true if you plan to use the external direct solvers such as SuperLU, MUMPS or Spooles. 2775 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2776 2777 Notes: 2778 If nnz is given then nz is ignored 2779 2780 The AIJ format (also called the Yale sparse matrix format or 2781 compressed row storage), is fully compatible with standard Fortran 77 2782 storage. That is, the stored row and column indices can begin at 2783 either one (as in Fortran) or zero. See the users' manual for details. 2784 2785 Specify the preallocated storage with either nz or nnz (not both). 2786 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2787 allocation. For large problems you MUST preallocate memory or you 2788 will get TERRIBLE performance, see the users' manual chapter on matrices. 2789 2790 By default, this format uses inodes (identical nodes) when possible, to 2791 improve numerical efficiency of matrix-vector products and solves. We 2792 search for consecutive rows with the same nonzero structure, thereby 2793 reusing matrix information to achieve increased efficiency. 2794 2795 Options Database Keys: 2796 + -mat_no_inode - Do not use inodes 2797 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2798 - -mat_aij_oneindex - Internally use indexing starting at 1 2799 rather than 0. Note that when calling MatSetValues(), 2800 the user still MUST index entries starting at 0! 2801 2802 Level: intermediate 2803 2804 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2805 2806 @*/ 2807 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2808 { 2809 PetscErrorCode ierr; 2810 2811 PetscFunctionBegin; 2812 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2813 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2814 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2815 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2816 PetscFunctionReturn(0); 2817 } 2818 2819 #undef __FUNCT__ 2820 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2821 /*@C 2822 MatSeqAIJSetPreallocation - For good matrix assembly performance 2823 the user should preallocate the matrix storage by setting the parameter nz 2824 (or the array nnz). By setting these parameters accurately, performance 2825 during matrix assembly can be increased by more than a factor of 50. 2826 2827 Collective on MPI_Comm 2828 2829 Input Parameters: 2830 + B - The matrix 2831 . nz - number of nonzeros per row (same for all rows) 2832 - nnz - array containing the number of nonzeros in the various rows 2833 (possibly different for each row) or PETSC_NULL 2834 2835 Notes: 2836 If nnz is given then nz is ignored 2837 2838 The AIJ format (also called the Yale sparse matrix format or 2839 compressed row storage), is fully compatible with standard Fortran 77 2840 storage. That is, the stored row and column indices can begin at 2841 either one (as in Fortran) or zero. See the users' manual for details. 2842 2843 Specify the preallocated storage with either nz or nnz (not both). 2844 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2845 allocation. For large problems you MUST preallocate memory or you 2846 will get TERRIBLE performance, see the users' manual chapter on matrices. 2847 2848 You can call MatGetInfo() to get information on how effective the preallocation was; 2849 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2850 You can also run with the option -info and look for messages with the string 2851 malloc in them to see if additional memory allocation was needed. 2852 2853 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2854 entries or columns indices 2855 2856 By default, this format uses inodes (identical nodes) when possible, to 2857 improve numerical efficiency of matrix-vector products and solves. We 2858 search for consecutive rows with the same nonzero structure, thereby 2859 reusing matrix information to achieve increased efficiency. 2860 2861 Options Database Keys: 2862 + -mat_no_inode - Do not use inodes 2863 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2864 - -mat_aij_oneindex - Internally use indexing starting at 1 2865 rather than 0. Note that when calling MatSetValues(), 2866 the user still MUST index entries starting at 0! 2867 2868 Level: intermediate 2869 2870 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2871 2872 @*/ 2873 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2874 { 2875 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2876 2877 PetscFunctionBegin; 2878 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2879 if (f) { 2880 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2881 } 2882 PetscFunctionReturn(0); 2883 } 2884 2885 EXTERN_C_BEGIN 2886 #undef __FUNCT__ 2887 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2888 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2889 { 2890 Mat_SeqAIJ *b; 2891 PetscTruth skipallocation = PETSC_FALSE; 2892 PetscErrorCode ierr; 2893 PetscInt i; 2894 2895 PetscFunctionBegin; 2896 2897 if (nz == MAT_SKIP_ALLOCATION) { 2898 skipallocation = PETSC_TRUE; 2899 nz = 0; 2900 } 2901 2902 B->rmap.bs = B->cmap.bs = 1; 2903 ierr = PetscMapSetUp(&B->rmap);CHKERRQ(ierr); 2904 ierr = PetscMapSetUp(&B->cmap);CHKERRQ(ierr); 2905 2906 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2907 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2908 if (nnz) { 2909 for (i=0; i<B->rmap.n; i++) { 2910 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2911 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); 2912 } 2913 } 2914 2915 B->preallocated = PETSC_TRUE; 2916 b = (Mat_SeqAIJ*)B->data; 2917 2918 if (!skipallocation) { 2919 if (!b->imax) { 2920 ierr = PetscMalloc2(B->rmap.n,PetscInt,&b->imax,B->rmap.n,PetscInt,&b->ilen);CHKERRQ(ierr); 2921 ierr = PetscLogObjectMemory(B,2*B->rmap.n*sizeof(PetscInt));CHKERRQ(ierr); 2922 } 2923 if (!nnz) { 2924 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2925 else if (nz <= 0) nz = 1; 2926 for (i=0; i<B->rmap.n; i++) b->imax[i] = nz; 2927 nz = nz*B->rmap.n; 2928 } else { 2929 nz = 0; 2930 for (i=0; i<B->rmap.n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2931 } 2932 /* b->ilen will count nonzeros in each row so far. */ 2933 for (i=0; i<B->rmap.n; i++) { b->ilen[i] = 0; } 2934 2935 /* allocate the matrix space */ 2936 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 2937 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap.n+1,PetscInt,&b->i);CHKERRQ(ierr); 2938 ierr = PetscLogObjectMemory(B,(B->rmap.n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 2939 b->i[0] = 0; 2940 for (i=1; i<B->rmap.n+1; i++) { 2941 b->i[i] = b->i[i-1] + b->imax[i-1]; 2942 } 2943 b->singlemalloc = PETSC_TRUE; 2944 b->free_a = PETSC_TRUE; 2945 b->free_ij = PETSC_TRUE; 2946 } else { 2947 b->free_a = PETSC_FALSE; 2948 b->free_ij = PETSC_FALSE; 2949 } 2950 2951 b->nz = 0; 2952 b->maxnz = nz; 2953 B->info.nz_unneeded = (double)b->maxnz; 2954 PetscFunctionReturn(0); 2955 } 2956 EXTERN_C_END 2957 2958 #undef __FUNCT__ 2959 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 2960 /*@ 2961 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 2962 2963 Input Parameters: 2964 + B - the matrix 2965 . i - the indices into j for the start of each row (starts with zero) 2966 . j - the column indices for each row (starts with zero) these must be sorted for each row 2967 - v - optional values in the matrix 2968 2969 Contributed by: Lisandro Dalchin 2970 2971 Level: developer 2972 2973 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 2974 2975 .keywords: matrix, aij, compressed row, sparse, sequential 2976 2977 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 2978 @*/ 2979 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 2980 { 2981 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 2982 PetscErrorCode ierr; 2983 2984 PetscFunctionBegin; 2985 PetscValidHeaderSpecific(B,MAT_COOKIE,1); 2986 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 2987 if (f) { 2988 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 2989 } 2990 PetscFunctionReturn(0); 2991 } 2992 2993 EXTERN_C_BEGIN 2994 #undef __FUNCT__ 2995 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 2996 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 2997 { 2998 PetscInt i; 2999 PetscInt m,n; 3000 PetscInt nz; 3001 PetscInt *nnz, nz_max = 0; 3002 PetscScalar *values; 3003 PetscErrorCode ierr; 3004 3005 PetscFunctionBegin; 3006 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3007 3008 if (Ii[0]) { 3009 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3010 } 3011 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3012 for(i = 0; i < m; i++) { 3013 nz = Ii[i+1]- Ii[i]; 3014 nz_max = PetscMax(nz_max, nz); 3015 if (nz < 0) { 3016 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3017 } 3018 nnz[i] = nz; 3019 } 3020 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3021 ierr = PetscFree(nnz);CHKERRQ(ierr); 3022 3023 if (v) { 3024 values = (PetscScalar*) v; 3025 } else { 3026 ierr = PetscMalloc((nz_max+1)*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3027 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3028 } 3029 3030 for(i = 0; i < m; i++) { 3031 nz = Ii[i+1] - Ii[i]; 3032 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3033 } 3034 3035 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3036 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3037 3038 if (!v) { 3039 ierr = PetscFree(values);CHKERRQ(ierr); 3040 } 3041 PetscFunctionReturn(0); 3042 } 3043 EXTERN_C_END 3044 3045 #include "src/mat/impls/dense/seq/dense.h" 3046 #include "src/inline/axpy.h" 3047 3048 #undef __FUNCT__ 3049 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3050 /* 3051 Computes (B'*A')' since computing B*A directly is untenable 3052 3053 n p p 3054 ( ) ( ) ( ) 3055 m ( A ) * n ( B ) = m ( C ) 3056 ( ) ( ) ( ) 3057 3058 */ 3059 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3060 { 3061 PetscErrorCode ierr; 3062 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3063 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3064 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3065 PetscInt i,n,m,q,p; 3066 const PetscInt *ii,*idx; 3067 const PetscScalar *b,*a,*a_q; 3068 PetscScalar *c,*c_q; 3069 3070 PetscFunctionBegin; 3071 m = A->rmap.n; 3072 n = A->cmap.n; 3073 p = B->cmap.n; 3074 a = sub_a->v; 3075 b = sub_b->a; 3076 c = sub_c->v; 3077 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3078 3079 ii = sub_b->i; 3080 idx = sub_b->j; 3081 for (i=0; i<n; i++) { 3082 q = ii[i+1] - ii[i]; 3083 while (q-->0) { 3084 c_q = c + m*(*idx); 3085 a_q = a + m*i; 3086 APXY(c_q,*b,a_q,m); 3087 idx++; 3088 b++; 3089 } 3090 } 3091 PetscFunctionReturn(0); 3092 } 3093 3094 #undef __FUNCT__ 3095 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3096 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3097 { 3098 PetscErrorCode ierr; 3099 PetscInt m=A->rmap.n,n=B->cmap.n; 3100 Mat Cmat; 3101 3102 PetscFunctionBegin; 3103 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); 3104 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3105 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3106 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3107 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3108 Cmat->assembled = PETSC_TRUE; 3109 *C = Cmat; 3110 PetscFunctionReturn(0); 3111 } 3112 3113 /* ----------------------------------------------------------------*/ 3114 #undef __FUNCT__ 3115 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3116 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3117 { 3118 PetscErrorCode ierr; 3119 3120 PetscFunctionBegin; 3121 if (scall == MAT_INITIAL_MATRIX){ 3122 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3123 } 3124 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3125 PetscFunctionReturn(0); 3126 } 3127 3128 3129 /*MC 3130 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3131 based on compressed sparse row format. 3132 3133 Options Database Keys: 3134 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3135 3136 Level: beginner 3137 3138 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3139 M*/ 3140 3141 EXTERN_C_BEGIN 3142 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3143 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3144 #if defined(PETSC_HAVE_MUMPS) 3145 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_mumps(Mat,MatFactorType,Mat*); 3146 #endif 3147 #if defined(PETSC_HAVE_SUPERLU) 3148 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3149 #endif 3150 #if defined(PETSC_HAVE_SPOOLES) 3151 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3152 #endif 3153 EXTERN_C_END 3154 3155 3156 EXTERN_C_BEGIN 3157 #undef __FUNCT__ 3158 #define __FUNCT__ "MatCreate_SeqAIJ" 3159 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3160 { 3161 Mat_SeqAIJ *b; 3162 PetscErrorCode ierr; 3163 PetscMPIInt size; 3164 3165 PetscFunctionBegin; 3166 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3167 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3168 3169 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3170 B->data = (void*)b; 3171 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3172 B->mapping = 0; 3173 b->row = 0; 3174 b->col = 0; 3175 b->icol = 0; 3176 b->reallocs = 0; 3177 b->ignorezeroentries = PETSC_FALSE; 3178 b->roworiented = PETSC_TRUE; 3179 b->nonew = 0; 3180 b->diag = 0; 3181 b->solve_work = 0; 3182 B->spptr = 0; 3183 b->saved_values = 0; 3184 b->idiag = 0; 3185 b->mdiag = 0; 3186 b->ssor_work = 0; 3187 b->omega = 1.0; 3188 b->fshift = 0.0; 3189 b->idiagvalid = PETSC_FALSE; 3190 b->keepzeroedrows = PETSC_FALSE; 3191 b->xtoy = 0; 3192 b->XtoY = 0; 3193 b->compressedrow.use = PETSC_FALSE; 3194 b->compressedrow.nrows = B->rmap.n; 3195 b->compressedrow.i = PETSC_NULL; 3196 b->compressedrow.rindex = PETSC_NULL; 3197 b->compressedrow.checked = PETSC_FALSE; 3198 B->same_nonzero = PETSC_FALSE; 3199 3200 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3201 #if defined(PETSC_HAVE_SUPERLU) 3202 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqaij_superlu_C", 3203 "MatGetFactor_seqaij_superlu", 3204 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3205 #endif 3206 #if defined(PETSC_HAVE_SPOOLES) 3207 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqaij_spooles_C", 3208 "MatGetFactor_seqaij_spooles", 3209 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3210 #endif 3211 #if defined(PETSC_HAVE_MUMPS) 3212 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqaij_mumps_C", 3213 "MatGetFactor_seqaij_mumps", 3214 MatGetFactor_seqaij_mumps);CHKERRQ(ierr); 3215 #endif 3216 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqaij_petsc_C", 3217 "MatGetFactor_seqaij_petsc", 3218 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3219 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3220 "MatSeqAIJSetColumnIndices_SeqAIJ", 3221 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3222 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3223 "MatStoreValues_SeqAIJ", 3224 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3225 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3226 "MatRetrieveValues_SeqAIJ", 3227 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3228 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3229 "MatConvert_SeqAIJ_SeqSBAIJ", 3230 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3231 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3232 "MatConvert_SeqAIJ_SeqBAIJ", 3233 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3234 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3235 "MatConvert_SeqAIJ_SeqCSRPERM", 3236 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3237 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3238 "MatConvert_SeqAIJ_SeqCRL", 3239 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3240 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3241 "MatIsTranspose_SeqAIJ", 3242 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3243 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3244 "MatIsHermitianTranspose_SeqAIJ", 3245 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3246 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3247 "MatSeqAIJSetPreallocation_SeqAIJ", 3248 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3249 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3250 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3251 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3252 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3253 "MatReorderForNonzeroDiagonal_SeqAIJ", 3254 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3255 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3256 "MatMatMult_SeqDense_SeqAIJ", 3257 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3258 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3259 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3260 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3261 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3262 "MatMatMultNumeric_SeqDense_SeqAIJ", 3263 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3264 ierr = MatCreate_Inode(B);CHKERRQ(ierr); 3265 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3266 PetscFunctionReturn(0); 3267 } 3268 EXTERN_C_END 3269 3270 #undef __FUNCT__ 3271 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3272 /* 3273 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3274 */ 3275 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3276 { 3277 Mat C = *B; 3278 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3279 PetscErrorCode ierr; 3280 PetscInt i,m = A->rmap.n; 3281 3282 PetscFunctionBegin; 3283 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3284 3285 c = (Mat_SeqAIJ*)C->data; 3286 3287 C->factor = A->factor; 3288 3289 c->row = 0; 3290 c->col = 0; 3291 c->icol = 0; 3292 c->reallocs = 0; 3293 3294 C->assembled = PETSC_TRUE; 3295 3296 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3297 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3298 for (i=0; i<m; i++) { 3299 c->imax[i] = a->imax[i]; 3300 c->ilen[i] = a->ilen[i]; 3301 } 3302 3303 /* allocate the matrix space */ 3304 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3305 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3306 c->singlemalloc = PETSC_TRUE; 3307 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3308 if (m > 0) { 3309 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3310 if (cpvalues == MAT_COPY_VALUES) { 3311 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3312 } else { 3313 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3314 } 3315 } 3316 3317 c->ignorezeroentries = a->ignorezeroentries; 3318 c->roworiented = a->roworiented; 3319 c->nonew = a->nonew; 3320 if (a->diag) { 3321 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3322 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3323 for (i=0; i<m; i++) { 3324 c->diag[i] = a->diag[i]; 3325 } 3326 } else c->diag = 0; 3327 c->solve_work = 0; 3328 c->saved_values = 0; 3329 c->idiag = 0; 3330 c->ssor_work = 0; 3331 c->keepzeroedrows = a->keepzeroedrows; 3332 c->free_a = PETSC_TRUE; 3333 c->free_ij = PETSC_TRUE; 3334 c->xtoy = 0; 3335 c->XtoY = 0; 3336 3337 c->nz = a->nz; 3338 c->maxnz = a->maxnz; 3339 C->preallocated = PETSC_TRUE; 3340 3341 c->compressedrow.use = a->compressedrow.use; 3342 c->compressedrow.nrows = a->compressedrow.nrows; 3343 c->compressedrow.checked = a->compressedrow.checked; 3344 if ( a->compressedrow.checked && a->compressedrow.use){ 3345 i = a->compressedrow.nrows; 3346 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 3347 c->compressedrow.rindex = c->compressedrow.i + i + 1; 3348 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3349 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3350 } else { 3351 c->compressedrow.use = PETSC_FALSE; 3352 c->compressedrow.i = PETSC_NULL; 3353 c->compressedrow.rindex = PETSC_NULL; 3354 } 3355 C->same_nonzero = A->same_nonzero; 3356 ierr = MatDuplicate_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3357 3358 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3359 PetscFunctionReturn(0); 3360 } 3361 3362 #undef __FUNCT__ 3363 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3364 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3365 { 3366 PetscErrorCode ierr; 3367 PetscInt n = A->rmap.n; 3368 3369 PetscFunctionBegin; 3370 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3371 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 3372 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3373 ierr = MatDuplicateNoCreate_SeqAIJ(A,cpvalues,B);CHKERRQ(ierr); 3374 PetscFunctionReturn(0); 3375 } 3376 3377 #undef __FUNCT__ 3378 #define __FUNCT__ "MatLoad_SeqAIJ" 3379 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, const MatType type,Mat *A) 3380 { 3381 Mat_SeqAIJ *a; 3382 Mat B; 3383 PetscErrorCode ierr; 3384 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3385 int fd; 3386 PetscMPIInt size; 3387 MPI_Comm comm; 3388 3389 PetscFunctionBegin; 3390 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3391 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3392 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 3393 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3394 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3395 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3396 M = header[1]; N = header[2]; nz = header[3]; 3397 3398 if (nz < 0) { 3399 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3400 } 3401 3402 /* read in row lengths */ 3403 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3404 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3405 3406 /* check if sum of rowlengths is same as nz */ 3407 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3408 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3409 3410 /* create our matrix */ 3411 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3412 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3413 ierr = MatSetType(B,type);CHKERRQ(ierr); 3414 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3415 a = (Mat_SeqAIJ*)B->data; 3416 3417 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3418 3419 /* read in nonzero values */ 3420 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3421 3422 /* set matrix "i" values */ 3423 a->i[0] = 0; 3424 for (i=1; i<= M; i++) { 3425 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3426 a->ilen[i-1] = rowlengths[i-1]; 3427 } 3428 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3429 3430 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3431 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3432 *A = B; 3433 PetscFunctionReturn(0); 3434 } 3435 3436 #undef __FUNCT__ 3437 #define __FUNCT__ "MatEqual_SeqAIJ" 3438 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3439 { 3440 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3441 PetscErrorCode ierr; 3442 3443 PetscFunctionBegin; 3444 /* If the matrix dimensions are not equal,or no of nonzeros */ 3445 if ((A->rmap.n != B->rmap.n) || (A->cmap.n != B->cmap.n) ||(a->nz != b->nz)) { 3446 *flg = PETSC_FALSE; 3447 PetscFunctionReturn(0); 3448 } 3449 3450 /* if the a->i are the same */ 3451 ierr = PetscMemcmp(a->i,b->i,(A->rmap.n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3452 if (!*flg) PetscFunctionReturn(0); 3453 3454 /* if a->j are the same */ 3455 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3456 if (!*flg) PetscFunctionReturn(0); 3457 3458 /* if a->a are the same */ 3459 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3460 3461 PetscFunctionReturn(0); 3462 3463 } 3464 3465 #undef __FUNCT__ 3466 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3467 /*@ 3468 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3469 provided by the user. 3470 3471 Collective on MPI_Comm 3472 3473 Input Parameters: 3474 + comm - must be an MPI communicator of size 1 3475 . m - number of rows 3476 . n - number of columns 3477 . i - row indices 3478 . j - column indices 3479 - a - matrix values 3480 3481 Output Parameter: 3482 . mat - the matrix 3483 3484 Level: intermediate 3485 3486 Notes: 3487 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3488 once the matrix is destroyed 3489 3490 You cannot set new nonzero locations into this matrix, that will generate an error. 3491 3492 The i and j indices are 0 based 3493 3494 The format which is used for the sparse matrix input, is equivalent to a 3495 row-major ordering.. i.e for the following matrix, the input data expected is 3496 as shown: 3497 3498 1 0 0 3499 2 0 3 3500 4 5 6 3501 3502 i = {0,1,3,6} [size = nrow+1 = 3+1] 3503 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3504 v = {1,2,3,4,5,6} [size = nz = 6] 3505 3506 3507 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3508 3509 @*/ 3510 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3511 { 3512 PetscErrorCode ierr; 3513 PetscInt ii; 3514 Mat_SeqAIJ *aij; 3515 #if defined(PETSC_USE_DEBUG) 3516 PetscInt jj; 3517 #endif 3518 3519 PetscFunctionBegin; 3520 if (i[0]) { 3521 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3522 } 3523 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3524 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3525 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3526 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3527 aij = (Mat_SeqAIJ*)(*mat)->data; 3528 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3529 3530 aij->i = i; 3531 aij->j = j; 3532 aij->a = a; 3533 aij->singlemalloc = PETSC_FALSE; 3534 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3535 aij->free_a = PETSC_FALSE; 3536 aij->free_ij = PETSC_FALSE; 3537 3538 for (ii=0; ii<m; ii++) { 3539 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3540 #if defined(PETSC_USE_DEBUG) 3541 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]); 3542 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3543 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); 3544 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); 3545 } 3546 #endif 3547 } 3548 #if defined(PETSC_USE_DEBUG) 3549 for (ii=0; ii<aij->i[m]; ii++) { 3550 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3551 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3552 } 3553 #endif 3554 3555 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3556 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3557 PetscFunctionReturn(0); 3558 } 3559 3560 #undef __FUNCT__ 3561 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3562 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3563 { 3564 PetscErrorCode ierr; 3565 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3566 3567 PetscFunctionBegin; 3568 if (coloring->ctype == IS_COLORING_GLOBAL) { 3569 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3570 a->coloring = coloring; 3571 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3572 PetscInt i,*larray; 3573 ISColoring ocoloring; 3574 ISColoringValue *colors; 3575 3576 /* set coloring for diagonal portion */ 3577 ierr = PetscMalloc((A->cmap.n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3578 for (i=0; i<A->cmap.n; i++) { 3579 larray[i] = i; 3580 } 3581 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap.n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3582 ierr = PetscMalloc((A->cmap.n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3583 for (i=0; i<A->cmap.n; i++) { 3584 colors[i] = coloring->colors[larray[i]]; 3585 } 3586 ierr = PetscFree(larray);CHKERRQ(ierr); 3587 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap.n,colors,&ocoloring);CHKERRQ(ierr); 3588 a->coloring = ocoloring; 3589 } 3590 PetscFunctionReturn(0); 3591 } 3592 3593 #if defined(PETSC_HAVE_ADIC) 3594 EXTERN_C_BEGIN 3595 #include "adic/ad_utils.h" 3596 EXTERN_C_END 3597 3598 #undef __FUNCT__ 3599 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3600 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3601 { 3602 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3603 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3604 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3605 ISColoringValue *color; 3606 3607 PetscFunctionBegin; 3608 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3609 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3610 color = a->coloring->colors; 3611 /* loop over rows */ 3612 for (i=0; i<m; i++) { 3613 nz = ii[i+1] - ii[i]; 3614 /* loop over columns putting computed value into matrix */ 3615 for (j=0; j<nz; j++) { 3616 *v++ = values[color[*jj++]]; 3617 } 3618 values += nlen; /* jump to next row of derivatives */ 3619 } 3620 PetscFunctionReturn(0); 3621 } 3622 #endif 3623 3624 #undef __FUNCT__ 3625 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3626 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3627 { 3628 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3629 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j; 3630 MatScalar *v = a->a; 3631 PetscScalar *values = (PetscScalar *)advalues; 3632 ISColoringValue *color; 3633 3634 PetscFunctionBegin; 3635 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3636 color = a->coloring->colors; 3637 /* loop over rows */ 3638 for (i=0; i<m; i++) { 3639 nz = ii[i+1] - ii[i]; 3640 /* loop over columns putting computed value into matrix */ 3641 for (j=0; j<nz; j++) { 3642 *v++ = values[color[*jj++]]; 3643 } 3644 values += nl; /* jump to next row of derivatives */ 3645 } 3646 PetscFunctionReturn(0); 3647 } 3648 3649 /* 3650 Special version for direct calls from Fortran 3651 */ 3652 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3653 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3654 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3655 #define matsetvaluesseqaij_ matsetvaluesseqaij 3656 #endif 3657 3658 /* Change these macros so can be used in void function */ 3659 #undef CHKERRQ 3660 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3661 #undef SETERRQ2 3662 #define SETERRQ2(ierr,b,c,d) CHKERRABORT(((PetscObject)A)->comm,ierr) 3663 3664 EXTERN_C_BEGIN 3665 #undef __FUNCT__ 3666 #define __FUNCT__ "matsetvaluesseqaij_" 3667 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3668 { 3669 Mat A = *AA; 3670 PetscInt m = *mm, n = *nn; 3671 InsertMode is = *isis; 3672 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3673 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3674 PetscInt *imax,*ai,*ailen; 3675 PetscErrorCode ierr; 3676 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3677 MatScalar *ap,value,*aa; 3678 PetscTruth ignorezeroentries = a->ignorezeroentries; 3679 PetscTruth roworiented = a->roworiented; 3680 3681 PetscFunctionBegin; 3682 MatPreallocated(A); 3683 imax = a->imax; 3684 ai = a->i; 3685 ailen = a->ilen; 3686 aj = a->j; 3687 aa = a->a; 3688 3689 for (k=0; k<m; k++) { /* loop over added rows */ 3690 row = im[k]; 3691 if (row < 0) continue; 3692 #if defined(PETSC_USE_DEBUG) 3693 if (row >= A->rmap.n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3694 #endif 3695 rp = aj + ai[row]; ap = aa + ai[row]; 3696 rmax = imax[row]; nrow = ailen[row]; 3697 low = 0; 3698 high = nrow; 3699 for (l=0; l<n; l++) { /* loop over added columns */ 3700 if (in[l] < 0) continue; 3701 #if defined(PETSC_USE_DEBUG) 3702 if (in[l] >= A->cmap.n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3703 #endif 3704 col = in[l]; 3705 if (roworiented) { 3706 value = v[l + k*n]; 3707 } else { 3708 value = v[k + l*m]; 3709 } 3710 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3711 3712 if (col <= lastcol) low = 0; else high = nrow; 3713 lastcol = col; 3714 while (high-low > 5) { 3715 t = (low+high)/2; 3716 if (rp[t] > col) high = t; 3717 else low = t; 3718 } 3719 for (i=low; i<high; i++) { 3720 if (rp[i] > col) break; 3721 if (rp[i] == col) { 3722 if (is == ADD_VALUES) ap[i] += value; 3723 else ap[i] = value; 3724 goto noinsert; 3725 } 3726 } 3727 if (value == 0.0 && ignorezeroentries) goto noinsert; 3728 if (nonew == 1) goto noinsert; 3729 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3730 MatSeqXAIJReallocateAIJ(A,A->rmap.n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3731 N = nrow++ - 1; a->nz++; high++; 3732 /* shift up all the later entries in this row */ 3733 for (ii=N; ii>=i; ii--) { 3734 rp[ii+1] = rp[ii]; 3735 ap[ii+1] = ap[ii]; 3736 } 3737 rp[i] = col; 3738 ap[i] = value; 3739 noinsert:; 3740 low = i + 1; 3741 } 3742 ailen[row] = nrow; 3743 } 3744 A->same_nonzero = PETSC_FALSE; 3745 PetscFunctionReturnVoid(); 3746 } 3747 EXTERN_C_END 3748