1 /*$Id: aij.c,v 1.385 2001/09/07 20:09:22 bsmith Exp $*/ 2 /* 3 Defines the basic matrix operations for the AIJ (compressed row) 4 matrix storage format. 5 */ 6 7 #include "src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 8 #include "src/vec/vecimpl.h" 9 #include "src/inline/spops.h" 10 #include "src/inline/dot.h" 11 #include "petscbt.h" 12 13 14 EXTERN int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 18 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int *ia[],int *ja[],PetscTruth *done) 19 { 20 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 21 int ierr,i,ishift; 22 23 PetscFunctionBegin; 24 *m = A->m; 25 if (!ia) PetscFunctionReturn(0); 26 ishift = 0; 27 if (symmetric && !A->structurally_symmetric) { 28 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 29 } else if (oshift == 1) { 30 int nz = a->i[A->m]; 31 /* malloc space and add 1 to i and j indices */ 32 ierr = PetscMalloc((A->m+1)*sizeof(int),ia);CHKERRQ(ierr); 33 ierr = PetscMalloc((nz+1)*sizeof(int),ja);CHKERRQ(ierr); 34 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 35 for (i=0; i<A->m+1; i++) (*ia)[i] = a->i[i] + 1; 36 } else { 37 *ia = a->i; *ja = a->j; 38 } 39 PetscFunctionReturn(0); 40 } 41 42 #undef __FUNCT__ 43 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 44 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int *ia[],int *ja[],PetscTruth *done) 45 { 46 int ierr; 47 48 PetscFunctionBegin; 49 if (!ia) PetscFunctionReturn(0); 50 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 51 ierr = PetscFree(*ia);CHKERRQ(ierr); 52 ierr = PetscFree(*ja);CHKERRQ(ierr); 53 } 54 PetscFunctionReturn(0); 55 } 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 59 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int *ia[],int *ja[],PetscTruth *done) 60 { 61 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 62 int ierr,i,*collengths,*cia,*cja,n = A->n,m = A->m; 63 int nz = a->i[m],row,*jj,mr,col; 64 65 PetscFunctionBegin; 66 *nn = A->n; 67 if (!ia) PetscFunctionReturn(0); 68 if (symmetric) { 69 ierr = MatToSymmetricIJ_SeqAIJ(A->m,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 70 } else { 71 ierr = PetscMalloc((n+1)*sizeof(int),&collengths);CHKERRQ(ierr); 72 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 73 ierr = PetscMalloc((n+1)*sizeof(int),&cia);CHKERRQ(ierr); 74 ierr = PetscMalloc((nz+1)*sizeof(int),&cja);CHKERRQ(ierr); 75 jj = a->j; 76 for (i=0; i<nz; i++) { 77 collengths[jj[i]]++; 78 } 79 cia[0] = oshift; 80 for (i=0; i<n; i++) { 81 cia[i+1] = cia[i] + collengths[i]; 82 } 83 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 84 jj = a->j; 85 for (row=0; row<m; row++) { 86 mr = a->i[row+1] - a->i[row]; 87 for (i=0; i<mr; i++) { 88 col = *jj++; 89 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 90 } 91 } 92 ierr = PetscFree(collengths);CHKERRQ(ierr); 93 *ia = cia; *ja = cja; 94 } 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 100 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int *ia[],int *ja[],PetscTruth *done) 101 { 102 int ierr; 103 104 PetscFunctionBegin; 105 if (!ia) PetscFunctionReturn(0); 106 107 ierr = PetscFree(*ia);CHKERRQ(ierr); 108 ierr = PetscFree(*ja);CHKERRQ(ierr); 109 110 PetscFunctionReturn(0); 111 } 112 113 #define CHUNKSIZE 15 114 115 #undef __FUNCT__ 116 #define __FUNCT__ "MatSetValues_SeqAIJ" 117 int MatSetValues_SeqAIJ(Mat A,int m,const int im[],int n,const int in[],const PetscScalar v[],InsertMode is) 118 { 119 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 120 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,sorted = a->sorted; 121 int *imax = a->imax,*ai = a->i,*ailen = a->ilen; 122 int *aj = a->j,nonew = a->nonew,ierr; 123 PetscScalar *ap,value,*aa = a->a; 124 PetscTruth ignorezeroentries = ((a->ignorezeroentries && is == ADD_VALUES) ? PETSC_TRUE:PETSC_FALSE); 125 PetscTruth roworiented = a->roworiented; 126 127 PetscFunctionBegin; 128 for (k=0; k<m; k++) { /* loop over added rows */ 129 row = im[k]; 130 if (row < 0) continue; 131 #if defined(PETSC_USE_BOPT_g) 132 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m-1); 133 #endif 134 rp = aj + ai[row]; ap = aa + ai[row]; 135 rmax = imax[row]; nrow = ailen[row]; 136 low = 0; 137 for (l=0; l<n; l++) { /* loop over added columns */ 138 if (in[l] < 0) continue; 139 #if defined(PETSC_USE_BOPT_g) 140 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->n-1); 141 #endif 142 col = in[l]; 143 if (roworiented) { 144 value = v[l + k*n]; 145 } else { 146 value = v[k + l*m]; 147 } 148 if (value == 0.0 && ignorezeroentries) continue; 149 150 if (!sorted) low = 0; high = nrow; 151 while (high-low > 5) { 152 t = (low+high)/2; 153 if (rp[t] > col) high = t; 154 else low = t; 155 } 156 for (i=low; i<high; i++) { 157 if (rp[i] > col) break; 158 if (rp[i] == col) { 159 if (is == ADD_VALUES) ap[i] += value; 160 else ap[i] = value; 161 goto noinsert; 162 } 163 } 164 if (nonew == 1) goto noinsert; 165 else if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%d,%d) in the matrix",row,col); 166 if (nrow >= rmax) { 167 /* there is no extra room in row, therefore enlarge */ 168 int new_nz = ai[A->m] + CHUNKSIZE,*new_i,*new_j; 169 size_t len; 170 PetscScalar *new_a; 171 172 if (nonew == -2) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%d,%d) in the matrix requiring new malloc()",row,col); 173 174 /* malloc new storage space */ 175 len = ((size_t) new_nz)*(sizeof(int)+sizeof(PetscScalar))+(A->m+1)*sizeof(int); 176 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 177 new_j = (int*)(new_a + new_nz); 178 new_i = new_j + new_nz; 179 180 /* copy over old data into new slots */ 181 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 182 for (ii=row+1; ii<A->m+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 183 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(int));CHKERRQ(ierr); 184 len = (((size_t) new_nz) - CHUNKSIZE - ai[row] - nrow ); 185 ierr = PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(int));CHKERRQ(ierr); 186 ierr = PetscMemcpy(new_a,aa,(((size_t) ai[row])+nrow)*sizeof(PetscScalar));CHKERRQ(ierr); 187 ierr = PetscMemcpy(new_a+ai[row]+nrow+CHUNKSIZE,aa+ai[row]+nrow,len*sizeof(PetscScalar));CHKERRQ(ierr); 188 /* free up old matrix storage */ 189 ierr = PetscFree(a->a);CHKERRQ(ierr); 190 if (!a->singlemalloc) { 191 ierr = PetscFree(a->i);CHKERRQ(ierr); 192 ierr = PetscFree(a->j);CHKERRQ(ierr); 193 } 194 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 195 a->singlemalloc = PETSC_TRUE; 196 197 rp = aj + ai[row]; ap = aa + ai[row] ; 198 rmax = imax[row] = imax[row] + CHUNKSIZE; 199 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(PetscScalar))); 200 a->maxnz += CHUNKSIZE; 201 a->reallocs++; 202 } 203 N = nrow++ - 1; a->nz++; 204 /* shift up all the later entries in this row */ 205 for (ii=N; ii>=i; ii--) { 206 rp[ii+1] = rp[ii]; 207 ap[ii+1] = ap[ii]; 208 } 209 rp[i] = col; 210 ap[i] = value; 211 noinsert:; 212 low = i + 1; 213 } 214 ailen[row] = nrow; 215 } 216 PetscFunctionReturn(0); 217 } 218 219 #undef __FUNCT__ 220 #define __FUNCT__ "MatGetValues_SeqAIJ" 221 int MatGetValues_SeqAIJ(Mat A,int m,const int im[],int n,const int in[],PetscScalar v[]) 222 { 223 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 224 int *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 225 int *ai = a->i,*ailen = a->ilen; 226 PetscScalar *ap,*aa = a->a,zero = 0.0; 227 228 PetscFunctionBegin; 229 for (k=0; k<m; k++) { /* loop over rows */ 230 row = im[k]; 231 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %d",row); 232 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %d max %d",row,A->m-1); 233 rp = aj + ai[row]; ap = aa + ai[row]; 234 nrow = ailen[row]; 235 for (l=0; l<n; l++) { /* loop over columns */ 236 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %d",in[l]); 237 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %d max %d",in[l],A->n-1); 238 col = in[l] ; 239 high = nrow; low = 0; /* assume unsorted */ 240 while (high-low > 5) { 241 t = (low+high)/2; 242 if (rp[t] > col) high = t; 243 else low = t; 244 } 245 for (i=low; i<high; i++) { 246 if (rp[i] > col) break; 247 if (rp[i] == col) { 248 *v++ = ap[i]; 249 goto finished; 250 } 251 } 252 *v++ = zero; 253 finished:; 254 } 255 } 256 PetscFunctionReturn(0); 257 } 258 259 260 #undef __FUNCT__ 261 #define __FUNCT__ "MatView_SeqAIJ_Binary" 262 int MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 263 { 264 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 265 int i,fd,*col_lens,ierr; 266 267 PetscFunctionBegin; 268 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 269 ierr = PetscMalloc((4+A->m)*sizeof(int),&col_lens);CHKERRQ(ierr); 270 col_lens[0] = MAT_FILE_COOKIE; 271 col_lens[1] = A->m; 272 col_lens[2] = A->n; 273 col_lens[3] = a->nz; 274 275 /* store lengths of each row and write (including header) to file */ 276 for (i=0; i<A->m; i++) { 277 col_lens[4+i] = a->i[i+1] - a->i[i]; 278 } 279 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,1);CHKERRQ(ierr); 280 ierr = PetscFree(col_lens);CHKERRQ(ierr); 281 282 /* store column indices (zero start index) */ 283 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,0);CHKERRQ(ierr); 284 285 /* store nonzero values */ 286 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,0);CHKERRQ(ierr); 287 PetscFunctionReturn(0); 288 } 289 290 extern int MatSeqAIJFactorInfo_UMFPACK(Mat,PetscViewer); 291 extern int MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 292 extern int MatFactorInfo_MUMPS(Mat,PetscViewer); 293 294 #undef __FUNCT__ 295 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 296 int MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 297 { 298 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 299 int ierr,i,j,m = A->m,shift=0; 300 char *name; 301 PetscViewerFormat format; 302 303 PetscFunctionBegin; 304 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 305 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 306 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_INFO) { 307 if (a->inode.size) { 308 ierr = PetscViewerASCIIPrintf(viewer,"using I-node routines: found %d nodes, limit used is %d\n",a->inode.node_count,a->inode.limit);CHKERRQ(ierr); 309 } else { 310 ierr = PetscViewerASCIIPrintf(viewer,"not using I-node routines\n");CHKERRQ(ierr); 311 } 312 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 313 int nofinalvalue = 0; 314 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->n-!shift)) { 315 nofinalvalue = 1; 316 } 317 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 318 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %d %d \n",m,A->n);CHKERRQ(ierr); 319 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %d \n",a->nz);CHKERRQ(ierr); 320 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%d,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 321 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 322 323 for (i=0; i<m; i++) { 324 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 325 #if defined(PETSC_USE_COMPLEX) 326 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); 327 #else 328 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 329 #endif 330 } 331 } 332 if (nofinalvalue) { 333 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e\n",m,A->n,0.0);CHKERRQ(ierr); 334 } 335 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 336 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 337 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 338 #if defined(PETSC_HAVE_UMFPACK) && !defined(PETSC_USE_SINGLE) && !defined(PETSC_USE_COMPLEX) 339 ierr = MatSeqAIJFactorInfo_UMFPACK(A,viewer);CHKERRQ(ierr); 340 #endif 341 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_SINGLE) && !defined(PETSC_USE_COMPLEX) 342 ierr = MatSeqAIJFactorInfo_Matlab(A,viewer);CHKERRQ(ierr); 343 #endif 344 #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_SINGLE) 345 ierr = MatFactorInfo_MUMPS(A,viewer);CHKERRQ(ierr); 346 #endif 347 PetscFunctionReturn(0); 348 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 349 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 350 for (i=0; i<m; i++) { 351 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 352 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 353 #if defined(PETSC_USE_COMPLEX) 354 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 355 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 356 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 357 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 358 } else if (PetscRealPart(a->a[j]) != 0.0) { 359 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 360 } 361 #else 362 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 363 #endif 364 } 365 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 366 } 367 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 368 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 369 int nzd=0,fshift=1,*sptr; 370 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 371 ierr = PetscMalloc((m+1)*sizeof(int),&sptr);CHKERRQ(ierr); 372 for (i=0; i<m; i++) { 373 sptr[i] = nzd+1; 374 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 375 if (a->j[j] >= i) { 376 #if defined(PETSC_USE_COMPLEX) 377 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 378 #else 379 if (a->a[j] != 0.0) nzd++; 380 #endif 381 } 382 } 383 } 384 sptr[m] = nzd+1; 385 ierr = PetscViewerASCIIPrintf(viewer," %d %d\n\n",m,nzd);CHKERRQ(ierr); 386 for (i=0; i<m+1; i+=6) { 387 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);} 388 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);} 389 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);} 390 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 391 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 392 else {ierr = PetscViewerASCIIPrintf(viewer," %d\n",sptr[i]);CHKERRQ(ierr);} 393 } 394 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 395 ierr = PetscFree(sptr);CHKERRQ(ierr); 396 for (i=0; i<m; i++) { 397 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 398 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %d ",a->j[j]+fshift);CHKERRQ(ierr);} 399 } 400 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 401 } 402 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 403 for (i=0; i<m; i++) { 404 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 405 if (a->j[j] >= i) { 406 #if defined(PETSC_USE_COMPLEX) 407 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 408 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 409 } 410 #else 411 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 412 #endif 413 } 414 } 415 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 416 } 417 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 418 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 419 int cnt = 0,jcnt; 420 PetscScalar value; 421 422 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 423 for (i=0; i<m; i++) { 424 jcnt = 0; 425 for (j=0; j<A->n; j++) { 426 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 427 value = a->a[cnt++]; 428 jcnt++; 429 } else { 430 value = 0.0; 431 } 432 #if defined(PETSC_USE_COMPLEX) 433 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 434 #else 435 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 436 #endif 437 } 438 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 439 } 440 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 441 } else { 442 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 443 for (i=0; i<m; i++) { 444 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 445 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 446 #if defined(PETSC_USE_COMPLEX) 447 if (PetscImaginaryPart(a->a[j]) > 0.0) { 448 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 449 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 450 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 451 } else { 452 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 453 } 454 #else 455 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 456 #endif 457 } 458 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 459 } 460 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 461 } 462 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 463 PetscFunctionReturn(0); 464 } 465 466 #undef __FUNCT__ 467 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 468 int MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 469 { 470 Mat A = (Mat) Aa; 471 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 472 int ierr,i,j,m = A->m,color; 473 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 474 PetscViewer viewer; 475 PetscViewerFormat format; 476 477 PetscFunctionBegin; 478 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 479 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 480 481 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 482 /* loop over matrix elements drawing boxes */ 483 484 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 485 /* Blue for negative, Cyan for zero and Red for positive */ 486 color = PETSC_DRAW_BLUE; 487 for (i=0; i<m; i++) { 488 y_l = m - i - 1.0; y_r = y_l + 1.0; 489 for (j=a->i[i]; j<a->i[i+1]; j++) { 490 x_l = a->j[j] ; x_r = x_l + 1.0; 491 #if defined(PETSC_USE_COMPLEX) 492 if (PetscRealPart(a->a[j]) >= 0.) continue; 493 #else 494 if (a->a[j] >= 0.) continue; 495 #endif 496 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 497 } 498 } 499 color = PETSC_DRAW_CYAN; 500 for (i=0; i<m; i++) { 501 y_l = m - i - 1.0; y_r = y_l + 1.0; 502 for (j=a->i[i]; j<a->i[i+1]; j++) { 503 x_l = a->j[j]; x_r = x_l + 1.0; 504 if (a->a[j] != 0.) continue; 505 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 506 } 507 } 508 color = PETSC_DRAW_RED; 509 for (i=0; i<m; i++) { 510 y_l = m - i - 1.0; y_r = y_l + 1.0; 511 for (j=a->i[i]; j<a->i[i+1]; j++) { 512 x_l = a->j[j]; x_r = x_l + 1.0; 513 #if defined(PETSC_USE_COMPLEX) 514 if (PetscRealPart(a->a[j]) <= 0.) continue; 515 #else 516 if (a->a[j] <= 0.) continue; 517 #endif 518 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 519 } 520 } 521 } else { 522 /* use contour shading to indicate magnitude of values */ 523 /* first determine max of all nonzero values */ 524 int nz = a->nz,count; 525 PetscDraw popup; 526 PetscReal scale; 527 528 for (i=0; i<nz; i++) { 529 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 530 } 531 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 532 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 533 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 534 count = 0; 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 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(a->a[count])); 540 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 541 count++; 542 } 543 } 544 } 545 PetscFunctionReturn(0); 546 } 547 548 #undef __FUNCT__ 549 #define __FUNCT__ "MatView_SeqAIJ_Draw" 550 int MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 551 { 552 int ierr; 553 PetscDraw draw; 554 PetscReal xr,yr,xl,yl,h,w; 555 PetscTruth isnull; 556 557 PetscFunctionBegin; 558 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 559 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 560 if (isnull) PetscFunctionReturn(0); 561 562 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 563 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 564 xr += w; yr += h; xl = -w; yl = -h; 565 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 566 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 567 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 568 PetscFunctionReturn(0); 569 } 570 571 #undef __FUNCT__ 572 #define __FUNCT__ "MatView_SeqAIJ" 573 int MatView_SeqAIJ(Mat A,PetscViewer viewer) 574 { 575 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 576 int ierr; 577 PetscTruth issocket,isascii,isbinary,isdraw; 578 579 PetscFunctionBegin; 580 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 581 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 582 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 583 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 584 if (issocket) { 585 ierr = PetscViewerSocketPutSparse_Private(viewer,A->m,A->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 586 } else if (isascii) { 587 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 588 } else if (isbinary) { 589 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 590 } else if (isdraw) { 591 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 592 } else { 593 SETERRQ1(1,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 594 } 595 PetscFunctionReturn(0); 596 } 597 598 EXTERN int Mat_AIJ_CheckInode(Mat,PetscTruth); 599 #undef __FUNCT__ 600 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 601 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 602 { 603 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 604 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax,ierr; 605 int m = A->m,*ip,N,*ailen = a->ilen,rmax = 0; 606 PetscScalar *aa = a->a,*ap; 607 #if defined(PETSC_HAVE_UMFPACK) || defined(PETSC_HAVE_MUMPS) 608 PetscTruth flag; 609 #endif 610 611 PetscFunctionBegin; 612 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 613 614 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 615 for (i=1; i<m; i++) { 616 /* move each row back by the amount of empty slots (fshift) before it*/ 617 fshift += imax[i-1] - ailen[i-1]; 618 rmax = PetscMax(rmax,ailen[i]); 619 if (fshift) { 620 ip = aj + ai[i] ; 621 ap = aa + ai[i] ; 622 N = ailen[i]; 623 for (j=0; j<N; j++) { 624 ip[j-fshift] = ip[j]; 625 ap[j-fshift] = ap[j]; 626 } 627 } 628 ai[i] = ai[i-1] + ailen[i-1]; 629 } 630 if (m) { 631 fshift += imax[m-1] - ailen[m-1]; 632 ai[m] = ai[m-1] + ailen[m-1]; 633 } 634 /* reset ilen and imax for each row */ 635 for (i=0; i<m; i++) { 636 ailen[i] = imax[i] = ai[i+1] - ai[i]; 637 } 638 a->nz = ai[m]; 639 640 /* diagonals may have moved, so kill the diagonal pointers */ 641 if (fshift && a->diag) { 642 ierr = PetscFree(a->diag);CHKERRQ(ierr); 643 PetscLogObjectMemory(A,-(m+1)*sizeof(int)); 644 a->diag = 0; 645 } 646 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded,%d used\n",m,A->n,fshift,a->nz); 647 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %d\n",a->reallocs); 648 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 649 a->reallocs = 0; 650 A->info.nz_unneeded = (double)fshift; 651 a->rmax = rmax; 652 653 /* check out for identical nodes. If found, use inode functions */ 654 ierr = Mat_AIJ_CheckInode(A,(PetscTruth)(!fshift));CHKERRQ(ierr); 655 656 #if defined(PETSC_HAVE_UMFPACK) 657 ierr = PetscOptionsHasName(A->prefix,"-mat_aij_umfpack",&flag);CHKERRQ(ierr); 658 if (flag) { ierr = MatUseUMFPACK_SeqAIJ(A);CHKERRQ(ierr); } 659 #endif 660 661 #if defined(PETSC_HAVE_MUMPS) 662 ierr = PetscOptionsHasName(A->prefix,"-mat_aij_mumps",&flag);CHKERRQ(ierr); 663 if (flag) { ierr = MatUseMUMPS_MPIAIJ(A);CHKERRQ(ierr); } 664 #endif 665 666 PetscFunctionReturn(0); 667 } 668 669 #undef __FUNCT__ 670 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 671 int MatZeroEntries_SeqAIJ(Mat A) 672 { 673 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 674 int ierr; 675 676 PetscFunctionBegin; 677 ierr = PetscMemzero(a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 678 PetscFunctionReturn(0); 679 } 680 681 #undef __FUNCT__ 682 #define __FUNCT__ "MatDestroy_SeqAIJ" 683 int MatDestroy_SeqAIJ(Mat A) 684 { 685 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 686 int ierr; 687 688 PetscFunctionBegin; 689 #if defined(PETSC_USE_LOG) 690 PetscLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",A->m,A->n,a->nz); 691 #endif 692 if (a->freedata) { 693 ierr = PetscFree(a->a);CHKERRQ(ierr); 694 if (!a->singlemalloc) { 695 ierr = PetscFree(a->i);CHKERRQ(ierr); 696 ierr = PetscFree(a->j);CHKERRQ(ierr); 697 } 698 } 699 if (a->row) { 700 ierr = ISDestroy(a->row);CHKERRQ(ierr); 701 } 702 if (a->col) { 703 ierr = ISDestroy(a->col);CHKERRQ(ierr); 704 } 705 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 706 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 707 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 708 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 709 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 710 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 711 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 712 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 713 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 714 if (a->xtoy) {ierr = PetscFree(a->xtoy);CHKERRQ(ierr);} 715 716 ierr = PetscFree(a);CHKERRQ(ierr); 717 PetscFunctionReturn(0); 718 } 719 720 #undef __FUNCT__ 721 #define __FUNCT__ "MatCompress_SeqAIJ" 722 int MatCompress_SeqAIJ(Mat A) 723 { 724 PetscFunctionBegin; 725 PetscFunctionReturn(0); 726 } 727 728 #undef __FUNCT__ 729 #define __FUNCT__ "MatSetOption_SeqAIJ" 730 int MatSetOption_SeqAIJ(Mat A,MatOption op) 731 { 732 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 733 734 PetscFunctionBegin; 735 switch (op) { 736 case MAT_ROW_ORIENTED: 737 a->roworiented = PETSC_TRUE; 738 break; 739 case MAT_KEEP_ZEROED_ROWS: 740 a->keepzeroedrows = PETSC_TRUE; 741 break; 742 case MAT_COLUMN_ORIENTED: 743 a->roworiented = PETSC_FALSE; 744 break; 745 case MAT_COLUMNS_SORTED: 746 a->sorted = PETSC_TRUE; 747 break; 748 case MAT_COLUMNS_UNSORTED: 749 a->sorted = PETSC_FALSE; 750 break; 751 case MAT_NO_NEW_NONZERO_LOCATIONS: 752 a->nonew = 1; 753 break; 754 case MAT_NEW_NONZERO_LOCATION_ERR: 755 a->nonew = -1; 756 break; 757 case MAT_NEW_NONZERO_ALLOCATION_ERR: 758 a->nonew = -2; 759 break; 760 case MAT_YES_NEW_NONZERO_LOCATIONS: 761 a->nonew = 0; 762 break; 763 case MAT_IGNORE_ZERO_ENTRIES: 764 a->ignorezeroentries = PETSC_TRUE; 765 break; 766 case MAT_USE_INODES: 767 a->inode.use = PETSC_TRUE; 768 break; 769 case MAT_DO_NOT_USE_INODES: 770 a->inode.use = PETSC_FALSE; 771 break; 772 case MAT_ROWS_SORTED: 773 case MAT_ROWS_UNSORTED: 774 case MAT_YES_NEW_DIAGONALS: 775 case MAT_IGNORE_OFF_PROC_ENTRIES: 776 case MAT_USE_HASH_TABLE: 777 PetscLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 778 break; 779 case MAT_NO_NEW_DIAGONALS: 780 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 781 case MAT_INODE_LIMIT_1: 782 a->inode.limit = 1; 783 break; 784 case MAT_INODE_LIMIT_2: 785 a->inode.limit = 2; 786 break; 787 case MAT_INODE_LIMIT_3: 788 a->inode.limit = 3; 789 break; 790 case MAT_INODE_LIMIT_4: 791 a->inode.limit = 4; 792 break; 793 case MAT_INODE_LIMIT_5: 794 a->inode.limit = 5; 795 break; 796 default: 797 SETERRQ(PETSC_ERR_SUP,"unknown option"); 798 } 799 PetscFunctionReturn(0); 800 } 801 802 #undef __FUNCT__ 803 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 804 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 805 { 806 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 807 int i,j,n,ierr; 808 PetscScalar *x,zero = 0.0; 809 810 PetscFunctionBegin; 811 ierr = VecSet(&zero,v);CHKERRQ(ierr); 812 ierr = VecGetArrayFast(v,&x);CHKERRQ(ierr); 813 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 814 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 815 for (i=0; i<A->m; i++) { 816 for (j=a->i[i]; j<a->i[i+1]; j++) { 817 if (a->j[j] == i) { 818 x[i] = a->a[j]; 819 break; 820 } 821 } 822 } 823 ierr = VecRestoreArrayFast(v,&x);CHKERRQ(ierr); 824 PetscFunctionReturn(0); 825 } 826 827 828 #undef __FUNCT__ 829 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 830 int MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 831 { 832 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 833 PetscScalar *x,*y; 834 int ierr,m = A->m; 835 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 836 PetscScalar *v,alpha; 837 int n,i,*idx; 838 #endif 839 840 PetscFunctionBegin; 841 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 842 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 843 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 844 845 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 846 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 847 #else 848 for (i=0; i<m; i++) { 849 idx = a->j + a->i[i] ; 850 v = a->a + a->i[i] ; 851 n = a->i[i+1] - a->i[i]; 852 alpha = x[i]; 853 while (n-->0) {y[*idx++] += alpha * *v++;} 854 } 855 #endif 856 PetscLogFlops(2*a->nz); 857 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 858 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 859 PetscFunctionReturn(0); 860 } 861 862 #undef __FUNCT__ 863 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 864 int MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 865 { 866 PetscScalar zero = 0.0; 867 int ierr; 868 869 PetscFunctionBegin; 870 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 871 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 872 PetscFunctionReturn(0); 873 } 874 875 876 #undef __FUNCT__ 877 #define __FUNCT__ "MatMult_SeqAIJ" 878 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 879 { 880 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 881 PetscScalar *x,*y,*v; 882 int ierr,m = A->m,*idx,*ii; 883 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 884 int n,i,jrow,j; 885 PetscScalar sum; 886 #endif 887 888 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 889 #pragma disjoint(*x,*y,*v) 890 #endif 891 892 PetscFunctionBegin; 893 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 894 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 895 idx = a->j; 896 v = a->a; 897 ii = a->i; 898 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 899 fortranmultaij_(&m,x,ii,idx,v,y); 900 #else 901 for (i=0; i<m; i++) { 902 jrow = ii[i]; 903 n = ii[i+1] - jrow; 904 sum = 0.0; 905 for (j=0; j<n; j++) { 906 sum += v[jrow]*x[idx[jrow]]; jrow++; 907 } 908 y[i] = sum; 909 } 910 #endif 911 PetscLogFlops(2*a->nz - m); 912 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 913 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 914 PetscFunctionReturn(0); 915 } 916 917 #undef __FUNCT__ 918 #define __FUNCT__ "MatMultAdd_SeqAIJ" 919 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 920 { 921 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 922 PetscScalar *x,*y,*z,*v; 923 int ierr,m = A->m,*idx,*ii; 924 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 925 int n,i,jrow,j; 926 PetscScalar sum; 927 #endif 928 929 PetscFunctionBegin; 930 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 931 ierr = VecGetArrayFast(yy,&y);CHKERRQ(ierr); 932 if (zz != yy) { 933 ierr = VecGetArrayFast(zz,&z);CHKERRQ(ierr); 934 } else { 935 z = y; 936 } 937 938 idx = a->j; 939 v = a->a; 940 ii = a->i; 941 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 942 fortranmultaddaij_(&m,x,ii,idx,v,y,z); 943 #else 944 for (i=0; i<m; i++) { 945 jrow = ii[i]; 946 n = ii[i+1] - jrow; 947 sum = y[i]; 948 for (j=0; j<n; j++) { 949 sum += v[jrow]*x[idx[jrow]]; jrow++; 950 } 951 z[i] = sum; 952 } 953 #endif 954 PetscLogFlops(2*a->nz); 955 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 956 ierr = VecRestoreArrayFast(yy,&y);CHKERRQ(ierr); 957 if (zz != yy) { 958 ierr = VecRestoreArrayFast(zz,&z);CHKERRQ(ierr); 959 } 960 PetscFunctionReturn(0); 961 } 962 963 /* 964 Adds diagonal pointers to sparse matrix structure. 965 */ 966 #undef __FUNCT__ 967 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 968 int MatMarkDiagonal_SeqAIJ(Mat A) 969 { 970 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 971 int i,j,*diag,m = A->m,ierr; 972 973 PetscFunctionBegin; 974 if (a->diag) PetscFunctionReturn(0); 975 976 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 977 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 978 for (i=0; i<A->m; i++) { 979 diag[i] = a->i[i+1]; 980 for (j=a->i[i]; j<a->i[i+1]; j++) { 981 if (a->j[j] == i) { 982 diag[i] = j; 983 break; 984 } 985 } 986 } 987 a->diag = diag; 988 PetscFunctionReturn(0); 989 } 990 991 /* 992 Checks for missing diagonals 993 */ 994 #undef __FUNCT__ 995 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 996 int MatMissingDiagonal_SeqAIJ(Mat A) 997 { 998 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 999 int *diag,*jj = a->j,i,ierr; 1000 1001 PetscFunctionBegin; 1002 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1003 diag = a->diag; 1004 for (i=0; i<A->m; i++) { 1005 if (jj[diag[i]] != i) { 1006 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 1007 } 1008 } 1009 PetscFunctionReturn(0); 1010 } 1011 1012 #undef __FUNCT__ 1013 #define __FUNCT__ "MatRelax_SeqAIJ" 1014 int MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,int lits,Vec xx) 1015 { 1016 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1017 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 1018 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 1019 int ierr,n = A->n,m = A->m,i; 1020 const int *idx,*diag; 1021 1022 PetscFunctionBegin; 1023 its = its*lits; 1024 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits); 1025 1026 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1027 diag = a->diag; 1028 if (!a->idiag) { 1029 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1030 a->ssor = a->idiag + m; 1031 mdiag = a->ssor + m; 1032 1033 v = a->a; 1034 1035 /* this is wrong when fshift omega changes each iteration */ 1036 if (omega == 1.0 && fshift == 0.0) { 1037 for (i=0; i<m; i++) { 1038 mdiag[i] = v[diag[i]]; 1039 a->idiag[i] = 1.0/v[diag[i]]; 1040 } 1041 PetscLogFlops(m); 1042 } else { 1043 for (i=0; i<m; i++) { 1044 mdiag[i] = v[diag[i]]; 1045 a->idiag[i] = omega/(fshift + v[diag[i]]); 1046 } 1047 PetscLogFlops(2*m); 1048 } 1049 } 1050 t = a->ssor; 1051 idiag = a->idiag; 1052 mdiag = a->idiag + 2*m; 1053 1054 ierr = VecGetArrayFast(xx,&x);CHKERRQ(ierr); 1055 if (xx != bb) { 1056 ierr = VecGetArrayFast(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1057 } else { 1058 b = x; 1059 } 1060 1061 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1062 xs = x; 1063 if (flag == SOR_APPLY_UPPER) { 1064 /* apply (U + D/omega) to the vector */ 1065 bs = b; 1066 for (i=0; i<m; i++) { 1067 d = fshift + a->a[diag[i]]; 1068 n = a->i[i+1] - diag[i] - 1; 1069 idx = a->j + diag[i] + 1; 1070 v = a->a + diag[i] + 1; 1071 sum = b[i]*d/omega; 1072 SPARSEDENSEDOT(sum,bs,v,idx,n); 1073 x[i] = sum; 1074 } 1075 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 1076 if (bb != xx) {ierr = VecRestoreArrayFast(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1077 PetscLogFlops(a->nz); 1078 PetscFunctionReturn(0); 1079 } 1080 1081 1082 /* Let A = L + U + D; where L is lower trianglar, 1083 U is upper triangular, E is diagonal; This routine applies 1084 1085 (L + E)^{-1} A (U + E)^{-1} 1086 1087 to a vector efficiently using Eisenstat's trick. This is for 1088 the case of SSOR preconditioner, so E is D/omega where omega 1089 is the relaxation factor. 1090 */ 1091 1092 if (flag == SOR_APPLY_LOWER) { 1093 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1094 } else if (flag & SOR_EISENSTAT) { 1095 /* Let A = L + U + D; where L is lower trianglar, 1096 U is upper triangular, E is diagonal; This routine applies 1097 1098 (L + E)^{-1} A (U + E)^{-1} 1099 1100 to a vector efficiently using Eisenstat's trick. This is for 1101 the case of SSOR preconditioner, so E is D/omega where omega 1102 is the relaxation factor. 1103 */ 1104 scale = (2.0/omega) - 1.0; 1105 1106 /* x = (E + U)^{-1} b */ 1107 for (i=m-1; i>=0; i--) { 1108 n = a->i[i+1] - diag[i] - 1; 1109 idx = a->j + diag[i] + 1; 1110 v = a->a + diag[i] + 1; 1111 sum = b[i]; 1112 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1113 x[i] = sum*idiag[i]; 1114 } 1115 1116 /* t = b - (2*E - D)x */ 1117 v = a->a; 1118 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1119 1120 /* t = (E + L)^{-1}t */ 1121 ts = t; 1122 diag = a->diag; 1123 for (i=0; i<m; i++) { 1124 n = diag[i] - a->i[i]; 1125 idx = a->j + a->i[i]; 1126 v = a->a + a->i[i]; 1127 sum = t[i]; 1128 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1129 t[i] = sum*idiag[i]; 1130 /* x = x + t */ 1131 x[i] += t[i]; 1132 } 1133 1134 PetscLogFlops(6*m-1 + 2*a->nz); 1135 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 1136 if (bb != xx) {ierr = VecRestoreArrayFast(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1137 PetscFunctionReturn(0); 1138 } 1139 if (flag & SOR_ZERO_INITIAL_GUESS) { 1140 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1141 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1142 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,idiag,a->a,(void*)b); 1143 #else 1144 for (i=0; i<m; i++) { 1145 n = diag[i] - a->i[i]; 1146 idx = a->j + a->i[i]; 1147 v = a->a + a->i[i]; 1148 sum = b[i]; 1149 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1150 x[i] = sum*idiag[i]; 1151 } 1152 #endif 1153 xb = x; 1154 PetscLogFlops(a->nz); 1155 } else xb = b; 1156 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1157 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1158 for (i=0; i<m; i++) { 1159 x[i] *= mdiag[i]; 1160 } 1161 PetscLogFlops(m); 1162 } 1163 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1164 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1165 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,idiag,a->a,(void*)xb); 1166 #else 1167 for (i=m-1; i>=0; i--) { 1168 n = a->i[i+1] - diag[i] - 1; 1169 idx = a->j + diag[i] + 1; 1170 v = a->a + diag[i] + 1; 1171 sum = xb[i]; 1172 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1173 x[i] = sum*idiag[i]; 1174 } 1175 #endif 1176 PetscLogFlops(a->nz); 1177 } 1178 its--; 1179 } 1180 while (its--) { 1181 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1182 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1183 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,diag,a->a,(void*)b); 1184 #else 1185 for (i=0; i<m; i++) { 1186 d = fshift + a->a[diag[i]]; 1187 n = a->i[i+1] - a->i[i]; 1188 idx = a->j + a->i[i]; 1189 v = a->a + a->i[i]; 1190 sum = b[i]; 1191 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1192 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1193 } 1194 #endif 1195 PetscLogFlops(a->nz); 1196 } 1197 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1198 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1199 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,diag,a->a,(void*)b); 1200 #else 1201 for (i=m-1; i>=0; i--) { 1202 d = fshift + a->a[diag[i]]; 1203 n = a->i[i+1] - a->i[i]; 1204 idx = a->j + a->i[i]; 1205 v = a->a + a->i[i]; 1206 sum = b[i]; 1207 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1208 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1209 } 1210 #endif 1211 PetscLogFlops(a->nz); 1212 } 1213 } 1214 ierr = VecRestoreArrayFast(xx,&x);CHKERRQ(ierr); 1215 if (bb != xx) {ierr = VecRestoreArrayFast(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1216 PetscFunctionReturn(0); 1217 } 1218 1219 #undef __FUNCT__ 1220 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1221 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1222 { 1223 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1224 1225 PetscFunctionBegin; 1226 info->rows_global = (double)A->m; 1227 info->columns_global = (double)A->n; 1228 info->rows_local = (double)A->m; 1229 info->columns_local = (double)A->n; 1230 info->block_size = 1.0; 1231 info->nz_allocated = (double)a->maxnz; 1232 info->nz_used = (double)a->nz; 1233 info->nz_unneeded = (double)(a->maxnz - a->nz); 1234 info->assemblies = (double)A->num_ass; 1235 info->mallocs = (double)a->reallocs; 1236 info->memory = A->mem; 1237 if (A->factor) { 1238 info->fill_ratio_given = A->info.fill_ratio_given; 1239 info->fill_ratio_needed = A->info.fill_ratio_needed; 1240 info->factor_mallocs = A->info.factor_mallocs; 1241 } else { 1242 info->fill_ratio_given = 0; 1243 info->fill_ratio_needed = 0; 1244 info->factor_mallocs = 0; 1245 } 1246 PetscFunctionReturn(0); 1247 } 1248 1249 EXTERN int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,MatFactorInfo*,Mat*); 1250 EXTERN int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1251 EXTERN int MatLUFactor_SeqAIJ(Mat,IS,IS,MatFactorInfo*); 1252 EXTERN int MatSolve_SeqAIJ(Mat,Vec,Vec); 1253 EXTERN int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1254 EXTERN int MatSolveTranspose_SeqAIJ(Mat,Vec,Vec); 1255 EXTERN int MatSolveTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1256 1257 #undef __FUNCT__ 1258 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1259 int MatZeroRows_SeqAIJ(Mat A,IS is,const PetscScalar *diag) 1260 { 1261 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1262 int i,ierr,N,*rows,m = A->m - 1; 1263 1264 PetscFunctionBegin; 1265 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1266 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1267 if (a->keepzeroedrows) { 1268 for (i=0; i<N; i++) { 1269 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1270 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1271 } 1272 if (diag) { 1273 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1274 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1275 for (i=0; i<N; i++) { 1276 a->a[a->diag[rows[i]]] = *diag; 1277 } 1278 } 1279 } else { 1280 if (diag) { 1281 for (i=0; i<N; i++) { 1282 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1283 if (a->ilen[rows[i]] > 0) { 1284 a->ilen[rows[i]] = 1; 1285 a->a[a->i[rows[i]]] = *diag; 1286 a->j[a->i[rows[i]]] = rows[i]; 1287 } else { /* in case row was completely empty */ 1288 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1289 } 1290 } 1291 } else { 1292 for (i=0; i<N; i++) { 1293 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1294 a->ilen[rows[i]] = 0; 1295 } 1296 } 1297 } 1298 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1299 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1300 PetscFunctionReturn(0); 1301 } 1302 1303 #undef __FUNCT__ 1304 #define __FUNCT__ "MatGetRow_SeqAIJ" 1305 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1306 { 1307 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1308 int *itmp; 1309 1310 PetscFunctionBegin; 1311 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %d out of range",row); 1312 1313 *nz = a->i[row+1] - a->i[row]; 1314 if (v) *v = a->a + a->i[row]; 1315 if (idx) { 1316 itmp = a->j + a->i[row]; 1317 if (*nz) { 1318 *idx = itmp; 1319 } 1320 else *idx = 0; 1321 } 1322 PetscFunctionReturn(0); 1323 } 1324 1325 /* remove this function? */ 1326 #undef __FUNCT__ 1327 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1328 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1329 { 1330 /* Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1331 int ierr; */ 1332 1333 PetscFunctionBegin; 1334 /* if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} */ 1335 PetscFunctionReturn(0); 1336 } 1337 1338 #undef __FUNCT__ 1339 #define __FUNCT__ "MatNorm_SeqAIJ" 1340 int MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1341 { 1342 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1343 PetscScalar *v = a->a; 1344 PetscReal sum = 0.0; 1345 int i,j,ierr; 1346 1347 PetscFunctionBegin; 1348 if (type == NORM_FROBENIUS) { 1349 for (i=0; i<a->nz; i++) { 1350 #if defined(PETSC_USE_COMPLEX) 1351 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1352 #else 1353 sum += (*v)*(*v); v++; 1354 #endif 1355 } 1356 *nrm = sqrt(sum); 1357 } else if (type == NORM_1) { 1358 PetscReal *tmp; 1359 int *jj = a->j; 1360 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1361 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1362 *nrm = 0.0; 1363 for (j=0; j<a->nz; j++) { 1364 tmp[*jj++] += PetscAbsScalar(*v); v++; 1365 } 1366 for (j=0; j<A->n; j++) { 1367 if (tmp[j] > *nrm) *nrm = tmp[j]; 1368 } 1369 ierr = PetscFree(tmp);CHKERRQ(ierr); 1370 } else if (type == NORM_INFINITY) { 1371 *nrm = 0.0; 1372 for (j=0; j<A->m; j++) { 1373 v = a->a + a->i[j]; 1374 sum = 0.0; 1375 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1376 sum += PetscAbsScalar(*v); v++; 1377 } 1378 if (sum > *nrm) *nrm = sum; 1379 } 1380 } else { 1381 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1382 } 1383 PetscFunctionReturn(0); 1384 } 1385 1386 #undef __FUNCT__ 1387 #define __FUNCT__ "MatTranspose_SeqAIJ" 1388 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1389 { 1390 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1391 Mat C; 1392 int i,ierr,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1393 PetscScalar *array = a->a; 1394 1395 PetscFunctionBegin; 1396 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1397 ierr = PetscMalloc((1+A->n)*sizeof(int),&col);CHKERRQ(ierr); 1398 ierr = PetscMemzero(col,(1+A->n)*sizeof(int));CHKERRQ(ierr); 1399 1400 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1401 ierr = MatCreateSeqAIJ(A->comm,A->n,m,0,col,&C);CHKERRQ(ierr); 1402 ierr = PetscFree(col);CHKERRQ(ierr); 1403 for (i=0; i<m; i++) { 1404 len = ai[i+1]-ai[i]; 1405 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1406 array += len; 1407 aj += len; 1408 } 1409 1410 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1411 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1412 1413 if (B) { 1414 *B = C; 1415 } else { 1416 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1417 } 1418 PetscFunctionReturn(0); 1419 } 1420 1421 EXTERN_C_BEGIN 1422 #undef __FUNCT__ 1423 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1424 int MatIsSymmetric_SeqAIJ(Mat A,Mat B,PetscTruth *f) 1425 { 1426 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1427 int *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1428 int ma,na,mb,nb, i,ierr; 1429 1430 PetscFunctionBegin; 1431 bij = (Mat_SeqAIJ *) B->data; 1432 1433 ierr = MatGetSize(A,&ma,&na); CHKERRQ(ierr); 1434 ierr = MatGetSize(B,&mb,&nb); CHKERRQ(ierr); 1435 if (ma!=nb || na!=mb) 1436 SETERRQ(1,"Incompatible A/B sizes for symmetry test"); 1437 aii = aij->i; bii = bij->i; 1438 adx = aij->j; bdx = bij->j; 1439 va = aij->a; vb = bij->a; 1440 ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr); 1441 ierr = PetscMalloc(mb*sizeof(int),&bptr); CHKERRQ(ierr); 1442 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1443 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1444 1445 *f = PETSC_TRUE; 1446 for (i=0; i<ma; i++) { 1447 /*printf("row %d spans %d--%d; we start @ %d\n", 1448 i,idx[ii[i]],idx[ii[i+1]-1],idx[aptr[i]]);*/ 1449 while (aptr[i]<aii[i+1]) { 1450 int idc,idr; PetscScalar vc,vr; 1451 /* column/row index/value */ 1452 idc = adx[aptr[i]]; idr = bdx[bptr[idc]]; 1453 vc = va[aptr[i]]; vr = vb[bptr[idc]]; 1454 /*printf("comparing %d: (%d,%d)=%e to (%d,%d)=%e\n", 1455 aptr[i],i,idc,vc,idc,idr,vr);*/ 1456 if (i!=idr || vc!=vr) { 1457 *f = PETSC_FALSE; goto done; 1458 } else { 1459 aptr[i]++; if (B || i!=idc) bptr[idc]++; 1460 } 1461 } 1462 } 1463 done: 1464 ierr = PetscFree(aptr); CHKERRQ(ierr); 1465 if (B) { 1466 ierr = PetscFree(bptr); CHKERRQ(ierr); 1467 } 1468 1469 PetscFunctionReturn(0); 1470 } 1471 EXTERN_C_END 1472 1473 #undef __FUNCT__ 1474 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1475 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1476 { 1477 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1478 PetscScalar *l,*r,x,*v; 1479 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1480 1481 PetscFunctionBegin; 1482 if (ll) { 1483 /* The local size is used so that VecMPI can be passed to this routine 1484 by MatDiagonalScale_MPIAIJ */ 1485 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1486 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1487 ierr = VecGetArrayFast(ll,&l);CHKERRQ(ierr); 1488 v = a->a; 1489 for (i=0; i<m; i++) { 1490 x = l[i]; 1491 M = a->i[i+1] - a->i[i]; 1492 for (j=0; j<M; j++) { (*v++) *= x;} 1493 } 1494 ierr = VecRestoreArrayFast(ll,&l);CHKERRQ(ierr); 1495 PetscLogFlops(nz); 1496 } 1497 if (rr) { 1498 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1499 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1500 ierr = VecGetArrayFast(rr,&r);CHKERRQ(ierr); 1501 v = a->a; jj = a->j; 1502 for (i=0; i<nz; i++) { 1503 (*v++) *= r[*jj++]; 1504 } 1505 ierr = VecRestoreArrayFast(rr,&r);CHKERRQ(ierr); 1506 PetscLogFlops(nz); 1507 } 1508 PetscFunctionReturn(0); 1509 } 1510 1511 #undef __FUNCT__ 1512 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1513 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1514 { 1515 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1516 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1517 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1518 int *irow,*icol,nrows,ncols; 1519 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1520 PetscScalar *a_new,*mat_a; 1521 Mat C; 1522 PetscTruth stride; 1523 1524 PetscFunctionBegin; 1525 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1526 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1527 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1528 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1529 1530 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1531 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1532 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1533 1534 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1535 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1536 if (stride && step == 1) { 1537 /* special case of contiguous rows */ 1538 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1539 starts = lens + nrows; 1540 /* loop over new rows determining lens and starting points */ 1541 for (i=0; i<nrows; i++) { 1542 kstart = ai[irow[i]]; 1543 kend = kstart + ailen[irow[i]]; 1544 for (k=kstart; k<kend; k++) { 1545 if (aj[k] >= first) { 1546 starts[i] = k; 1547 break; 1548 } 1549 } 1550 sum = 0; 1551 while (k < kend) { 1552 if (aj[k++] >= first+ncols) break; 1553 sum++; 1554 } 1555 lens[i] = sum; 1556 } 1557 /* create submatrix */ 1558 if (scall == MAT_REUSE_MATRIX) { 1559 int n_cols,n_rows; 1560 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1561 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1562 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1563 C = *B; 1564 } else { 1565 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1566 } 1567 c = (Mat_SeqAIJ*)C->data; 1568 1569 /* loop over rows inserting into submatrix */ 1570 a_new = c->a; 1571 j_new = c->j; 1572 i_new = c->i; 1573 1574 for (i=0; i<nrows; i++) { 1575 ii = starts[i]; 1576 lensi = lens[i]; 1577 for (k=0; k<lensi; k++) { 1578 *j_new++ = aj[ii+k] - first; 1579 } 1580 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1581 a_new += lensi; 1582 i_new[i+1] = i_new[i] + lensi; 1583 c->ilen[i] = lensi; 1584 } 1585 ierr = PetscFree(lens);CHKERRQ(ierr); 1586 } else { 1587 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1588 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1589 1590 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1591 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1592 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1593 /* determine lens of each row */ 1594 for (i=0; i<nrows; i++) { 1595 kstart = ai[irow[i]]; 1596 kend = kstart + a->ilen[irow[i]]; 1597 lens[i] = 0; 1598 for (k=kstart; k<kend; k++) { 1599 if (smap[aj[k]]) { 1600 lens[i]++; 1601 } 1602 } 1603 } 1604 /* Create and fill new matrix */ 1605 if (scall == MAT_REUSE_MATRIX) { 1606 PetscTruth equal; 1607 1608 c = (Mat_SeqAIJ *)((*B)->data); 1609 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1610 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1611 if (!equal) { 1612 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1613 } 1614 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1615 C = *B; 1616 } else { 1617 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1618 } 1619 c = (Mat_SeqAIJ *)(C->data); 1620 for (i=0; i<nrows; i++) { 1621 row = irow[i]; 1622 kstart = ai[row]; 1623 kend = kstart + a->ilen[row]; 1624 mat_i = c->i[i]; 1625 mat_j = c->j + mat_i; 1626 mat_a = c->a + mat_i; 1627 mat_ilen = c->ilen + i; 1628 for (k=kstart; k<kend; k++) { 1629 if ((tcol=smap[a->j[k]])) { 1630 *mat_j++ = tcol - 1; 1631 *mat_a++ = a->a[k]; 1632 (*mat_ilen)++; 1633 1634 } 1635 } 1636 } 1637 /* Free work space */ 1638 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1639 ierr = PetscFree(smap);CHKERRQ(ierr); 1640 ierr = PetscFree(lens);CHKERRQ(ierr); 1641 } 1642 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1643 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1644 1645 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1646 *B = C; 1647 PetscFunctionReturn(0); 1648 } 1649 1650 /* 1651 */ 1652 #undef __FUNCT__ 1653 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1654 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1655 { 1656 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1657 int ierr; 1658 Mat outA; 1659 PetscTruth row_identity,col_identity; 1660 1661 PetscFunctionBegin; 1662 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1663 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1664 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1665 if (!row_identity || !col_identity) { 1666 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1667 } 1668 1669 outA = inA; 1670 inA->factor = FACTOR_LU; 1671 a->row = row; 1672 a->col = col; 1673 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1674 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1675 1676 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1677 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1678 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1679 PetscLogObjectParent(inA,a->icol); 1680 1681 if (!a->solve_work) { /* this matrix may have been factored before */ 1682 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1683 } 1684 1685 if (!a->diag) { 1686 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1687 } 1688 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1689 PetscFunctionReturn(0); 1690 } 1691 1692 #include "petscblaslapack.h" 1693 #undef __FUNCT__ 1694 #define __FUNCT__ "MatScale_SeqAIJ" 1695 int MatScale_SeqAIJ(const PetscScalar *alpha,Mat inA) 1696 { 1697 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1698 int one = 1; 1699 1700 PetscFunctionBegin; 1701 BLscal_(&a->nz,(PetscScalar*)alpha,a->a,&one); 1702 PetscLogFlops(a->nz); 1703 PetscFunctionReturn(0); 1704 } 1705 1706 #undef __FUNCT__ 1707 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1708 int MatGetSubMatrices_SeqAIJ(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1709 { 1710 int ierr,i; 1711 1712 PetscFunctionBegin; 1713 if (scall == MAT_INITIAL_MATRIX) { 1714 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1715 } 1716 1717 for (i=0; i<n; i++) { 1718 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1719 } 1720 PetscFunctionReturn(0); 1721 } 1722 1723 #undef __FUNCT__ 1724 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1725 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1726 { 1727 PetscFunctionBegin; 1728 *bs = 1; 1729 PetscFunctionReturn(0); 1730 } 1731 1732 #undef __FUNCT__ 1733 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1734 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS is[],int ov) 1735 { 1736 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1737 int row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1738 int start,end,*ai,*aj; 1739 PetscBT table; 1740 1741 PetscFunctionBegin; 1742 m = A->m; 1743 ai = a->i; 1744 aj = a->j; 1745 1746 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1747 1748 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1749 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1750 1751 for (i=0; i<is_max; i++) { 1752 /* Initialize the two local arrays */ 1753 isz = 0; 1754 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1755 1756 /* Extract the indices, assume there can be duplicate entries */ 1757 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1758 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1759 1760 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1761 for (j=0; j<n ; ++j){ 1762 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1763 } 1764 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1765 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1766 1767 k = 0; 1768 for (j=0; j<ov; j++){ /* for each overlap */ 1769 n = isz; 1770 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1771 row = nidx[k]; 1772 start = ai[row]; 1773 end = ai[row+1]; 1774 for (l = start; l<end ; l++){ 1775 val = aj[l] ; 1776 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1777 } 1778 } 1779 } 1780 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1781 } 1782 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1783 ierr = PetscFree(nidx);CHKERRQ(ierr); 1784 PetscFunctionReturn(0); 1785 } 1786 1787 /* -------------------------------------------------------------- */ 1788 #undef __FUNCT__ 1789 #define __FUNCT__ "MatPermute_SeqAIJ" 1790 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1791 { 1792 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1793 PetscScalar *vwork; 1794 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1795 int *row,*col,*cnew,j,*lens; 1796 IS icolp,irowp; 1797 1798 PetscFunctionBegin; 1799 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1800 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1801 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1802 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1803 1804 /* determine lengths of permuted rows */ 1805 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1806 for (i=0; i<m; i++) { 1807 lens[row[i]] = a->i[i+1] - a->i[i]; 1808 } 1809 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1810 ierr = PetscFree(lens);CHKERRQ(ierr); 1811 1812 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1813 for (i=0; i<m; i++) { 1814 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1815 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1816 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1817 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1818 } 1819 ierr = PetscFree(cnew);CHKERRQ(ierr); 1820 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1821 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1822 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1823 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1824 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1825 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1826 PetscFunctionReturn(0); 1827 } 1828 1829 #undef __FUNCT__ 1830 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1831 int MatPrintHelp_SeqAIJ(Mat A) 1832 { 1833 static PetscTruth called = PETSC_FALSE; 1834 MPI_Comm comm = A->comm; 1835 int ierr; 1836 1837 PetscFunctionBegin; 1838 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1839 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1840 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1841 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1842 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1843 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1844 #if defined(PETSC_HAVE_ESSL) 1845 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1846 #endif 1847 #if defined(PETSC_HAVE_LUSOL) 1848 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1849 #endif 1850 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1851 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1852 #endif 1853 PetscFunctionReturn(0); 1854 } 1855 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1856 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1857 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatFactorInfo*,IS,IS,Mat*); 1858 #undef __FUNCT__ 1859 #define __FUNCT__ "MatCopy_SeqAIJ" 1860 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1861 { 1862 int ierr; 1863 1864 PetscFunctionBegin; 1865 /* If the two matrices have the same copy implementation, use fast copy. */ 1866 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1867 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1868 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1869 1870 if (a->i[A->m] != b->i[B->m]) { 1871 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1872 } 1873 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1874 } else { 1875 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1876 } 1877 PetscFunctionReturn(0); 1878 } 1879 1880 #undef __FUNCT__ 1881 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1882 int MatSetUpPreallocation_SeqAIJ(Mat A) 1883 { 1884 int ierr; 1885 1886 PetscFunctionBegin; 1887 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1888 PetscFunctionReturn(0); 1889 } 1890 1891 #undef __FUNCT__ 1892 #define __FUNCT__ "MatGetArray_SeqAIJ" 1893 int MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1894 { 1895 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1896 PetscFunctionBegin; 1897 *array = a->a; 1898 PetscFunctionReturn(0); 1899 } 1900 1901 #undef __FUNCT__ 1902 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1903 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1904 { 1905 PetscFunctionBegin; 1906 PetscFunctionReturn(0); 1907 } 1908 1909 #undef __FUNCT__ 1910 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1911 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1912 { 1913 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1914 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1915 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1916 PetscScalar *vscale_array; 1917 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1918 Vec w1,w2,w3; 1919 void *fctx = coloring->fctx; 1920 PetscTruth flg; 1921 1922 PetscFunctionBegin; 1923 if (!coloring->w1) { 1924 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1925 PetscLogObjectParent(coloring,coloring->w1); 1926 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1927 PetscLogObjectParent(coloring,coloring->w2); 1928 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1929 PetscLogObjectParent(coloring,coloring->w3); 1930 } 1931 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1932 1933 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1934 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1935 if (flg) { 1936 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1937 } else { 1938 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1939 } 1940 1941 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1942 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1943 1944 /* 1945 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1946 coloring->F for the coarser grids from the finest 1947 */ 1948 if (coloring->F) { 1949 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1950 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1951 if (m1 != m2) { 1952 coloring->F = 0; 1953 } 1954 } 1955 1956 if (coloring->F) { 1957 w1 = coloring->F; 1958 coloring->F = 0; 1959 } else { 1960 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1961 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1962 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1963 } 1964 1965 /* 1966 Compute all the scale factors and share with other processors 1967 */ 1968 ierr = VecGetArrayFast(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1969 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1970 for (k=0; k<coloring->ncolors; k++) { 1971 /* 1972 Loop over each column associated with color adding the 1973 perturbation to the vector w3. 1974 */ 1975 for (l=0; l<coloring->ncolumns[k]; l++) { 1976 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1977 dx = xx[col]; 1978 if (dx == 0.0) dx = 1.0; 1979 #if !defined(PETSC_USE_COMPLEX) 1980 if (dx < umin && dx >= 0.0) dx = umin; 1981 else if (dx < 0.0 && dx > -umin) dx = -umin; 1982 #else 1983 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1984 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1985 #endif 1986 dx *= epsilon; 1987 vscale_array[col] = 1.0/dx; 1988 } 1989 } 1990 vscale_array = vscale_array + start;ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1991 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1992 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1993 1994 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1995 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 1996 1997 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 1998 else vscaleforrow = coloring->columnsforrow; 1999 2000 ierr = VecGetArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2001 /* 2002 Loop over each color 2003 */ 2004 for (k=0; k<coloring->ncolors; k++) { 2005 coloring->currentcolor = k; 2006 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2007 ierr = VecGetArrayFast(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2008 /* 2009 Loop over each column associated with color adding the 2010 perturbation to the vector w3. 2011 */ 2012 for (l=0; l<coloring->ncolumns[k]; l++) { 2013 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2014 dx = xx[col]; 2015 if (dx == 0.0) dx = 1.0; 2016 #if !defined(PETSC_USE_COMPLEX) 2017 if (dx < umin && dx >= 0.0) dx = umin; 2018 else if (dx < 0.0 && dx > -umin) dx = -umin; 2019 #else 2020 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2021 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2022 #endif 2023 dx *= epsilon; 2024 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2025 w3_array[col] += dx; 2026 } 2027 w3_array = w3_array + start; ierr = VecRestoreArrayFast(w3,&w3_array);CHKERRQ(ierr); 2028 2029 /* 2030 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2031 */ 2032 2033 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2034 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2035 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2036 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2037 2038 /* 2039 Loop over rows of vector, putting results into Jacobian matrix 2040 */ 2041 ierr = VecGetArrayFast(w2,&y);CHKERRQ(ierr); 2042 for (l=0; l<coloring->nrows[k]; l++) { 2043 row = coloring->rows[k][l]; 2044 col = coloring->columnsforrow[k][l]; 2045 y[row] *= vscale_array[vscaleforrow[k][l]]; 2046 srow = row + start; 2047 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2048 } 2049 ierr = VecRestoreArrayFast(w2,&y);CHKERRQ(ierr); 2050 } 2051 coloring->currentcolor = k; 2052 ierr = VecRestoreArrayFast(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2053 xx = xx + start; ierr = VecRestoreArrayFast(x1,&xx);CHKERRQ(ierr); 2054 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2055 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2056 PetscFunctionReturn(0); 2057 } 2058 2059 #include "petscblaslapack.h" 2060 #undef __FUNCT__ 2061 #define __FUNCT__ "MatAXPY_SeqAIJ" 2062 int MatAXPY_SeqAIJ(const PetscScalar *a,Mat X,Mat Y,MatStructure str) 2063 { 2064 int ierr,one=1,i; 2065 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2066 2067 PetscFunctionBegin; 2068 if (str == SAME_NONZERO_PATTERN) { 2069 BLaxpy_(&x->nz,(PetscScalar*)a,x->a,&one,y->a,&one); 2070 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2071 if (y->xtoy && y->XtoY != X) { 2072 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2073 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2074 } 2075 if (!y->xtoy) { /* get xtoy */ 2076 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2077 y->XtoY = X; 2078 } 2079 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2080 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2081 } else { 2082 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2083 } 2084 PetscFunctionReturn(0); 2085 } 2086 2087 /* -------------------------------------------------------------------*/ 2088 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2089 MatGetRow_SeqAIJ, 2090 MatRestoreRow_SeqAIJ, 2091 MatMult_SeqAIJ, 2092 MatMultAdd_SeqAIJ, 2093 MatMultTranspose_SeqAIJ, 2094 MatMultTransposeAdd_SeqAIJ, 2095 MatSolve_SeqAIJ, 2096 MatSolveAdd_SeqAIJ, 2097 MatSolveTranspose_SeqAIJ, 2098 MatSolveTransposeAdd_SeqAIJ, 2099 MatLUFactor_SeqAIJ, 2100 0, 2101 MatRelax_SeqAIJ, 2102 MatTranspose_SeqAIJ, 2103 MatGetInfo_SeqAIJ, 2104 MatEqual_SeqAIJ, 2105 MatGetDiagonal_SeqAIJ, 2106 MatDiagonalScale_SeqAIJ, 2107 MatNorm_SeqAIJ, 2108 0, 2109 MatAssemblyEnd_SeqAIJ, 2110 MatCompress_SeqAIJ, 2111 MatSetOption_SeqAIJ, 2112 MatZeroEntries_SeqAIJ, 2113 MatZeroRows_SeqAIJ, 2114 MatLUFactorSymbolic_SeqAIJ, 2115 MatLUFactorNumeric_SeqAIJ, 2116 MatCholeskyFactorSymbolic_SeqAIJ, 2117 MatCholeskyFactorNumeric_SeqAIJ, 2118 MatSetUpPreallocation_SeqAIJ, 2119 MatILUFactorSymbolic_SeqAIJ, 2120 MatICCFactorSymbolic_SeqAIJ, 2121 MatGetArray_SeqAIJ, 2122 MatRestoreArray_SeqAIJ, 2123 MatDuplicate_SeqAIJ, 2124 0, 2125 0, 2126 MatILUFactor_SeqAIJ, 2127 0, 2128 MatAXPY_SeqAIJ, 2129 MatGetSubMatrices_SeqAIJ, 2130 MatIncreaseOverlap_SeqAIJ, 2131 MatGetValues_SeqAIJ, 2132 MatCopy_SeqAIJ, 2133 MatPrintHelp_SeqAIJ, 2134 MatScale_SeqAIJ, 2135 0, 2136 0, 2137 MatILUDTFactor_SeqAIJ, 2138 MatGetBlockSize_SeqAIJ, 2139 MatGetRowIJ_SeqAIJ, 2140 MatRestoreRowIJ_SeqAIJ, 2141 MatGetColumnIJ_SeqAIJ, 2142 MatRestoreColumnIJ_SeqAIJ, 2143 MatFDColoringCreate_SeqAIJ, 2144 0, 2145 0, 2146 MatPermute_SeqAIJ, 2147 0, 2148 0, 2149 MatDestroy_SeqAIJ, 2150 MatView_SeqAIJ, 2151 MatGetPetscMaps_Petsc, 2152 0, 2153 0, 2154 0, 2155 0, 2156 0, 2157 0, 2158 0, 2159 0, 2160 MatSetColoring_SeqAIJ, 2161 MatSetValuesAdic_SeqAIJ, 2162 MatSetValuesAdifor_SeqAIJ, 2163 MatFDColoringApply_SeqAIJ}; 2164 2165 EXTERN_C_BEGIN 2166 #undef __FUNCT__ 2167 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2168 2169 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2170 { 2171 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2172 int i,nz,n; 2173 2174 PetscFunctionBegin; 2175 2176 nz = aij->maxnz; 2177 n = mat->n; 2178 for (i=0; i<nz; i++) { 2179 aij->j[i] = indices[i]; 2180 } 2181 aij->nz = nz; 2182 for (i=0; i<n; i++) { 2183 aij->ilen[i] = aij->imax[i]; 2184 } 2185 2186 PetscFunctionReturn(0); 2187 } 2188 EXTERN_C_END 2189 2190 #undef __FUNCT__ 2191 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2192 /*@ 2193 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2194 in the matrix. 2195 2196 Input Parameters: 2197 + mat - the SeqAIJ matrix 2198 - indices - the column indices 2199 2200 Level: advanced 2201 2202 Notes: 2203 This can be called if you have precomputed the nonzero structure of the 2204 matrix and want to provide it to the matrix object to improve the performance 2205 of the MatSetValues() operation. 2206 2207 You MUST have set the correct numbers of nonzeros per row in the call to 2208 MatCreateSeqAIJ(). 2209 2210 MUST be called before any calls to MatSetValues(); 2211 2212 The indices should start with zero, not one. 2213 2214 @*/ 2215 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2216 { 2217 int ierr,(*f)(Mat,int *); 2218 2219 PetscFunctionBegin; 2220 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2221 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2222 if (f) { 2223 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2224 } else { 2225 SETERRQ(1,"Wrong type of matrix to set column indices"); 2226 } 2227 PetscFunctionReturn(0); 2228 } 2229 2230 /* ----------------------------------------------------------------------------------------*/ 2231 2232 EXTERN_C_BEGIN 2233 #undef __FUNCT__ 2234 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2235 int MatStoreValues_SeqAIJ(Mat mat) 2236 { 2237 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2238 size_t nz = aij->i[mat->m],ierr; 2239 2240 PetscFunctionBegin; 2241 if (aij->nonew != 1) { 2242 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2243 } 2244 2245 /* allocate space for values if not already there */ 2246 if (!aij->saved_values) { 2247 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2248 } 2249 2250 /* copy values over */ 2251 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2252 PetscFunctionReturn(0); 2253 } 2254 EXTERN_C_END 2255 2256 #undef __FUNCT__ 2257 #define __FUNCT__ "MatStoreValues" 2258 /*@ 2259 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2260 example, reuse of the linear part of a Jacobian, while recomputing the 2261 nonlinear portion. 2262 2263 Collect on Mat 2264 2265 Input Parameters: 2266 . mat - the matrix (currently on AIJ matrices support this option) 2267 2268 Level: advanced 2269 2270 Common Usage, with SNESSolve(): 2271 $ Create Jacobian matrix 2272 $ Set linear terms into matrix 2273 $ Apply boundary conditions to matrix, at this time matrix must have 2274 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2275 $ boundary conditions again will not change the nonzero structure 2276 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2277 $ ierr = MatStoreValues(mat); 2278 $ Call SNESSetJacobian() with matrix 2279 $ In your Jacobian routine 2280 $ ierr = MatRetrieveValues(mat); 2281 $ Set nonlinear terms in matrix 2282 2283 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2284 $ // build linear portion of Jacobian 2285 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2286 $ ierr = MatStoreValues(mat); 2287 $ loop over nonlinear iterations 2288 $ ierr = MatRetrieveValues(mat); 2289 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2290 $ // call MatAssemblyBegin/End() on matrix 2291 $ Solve linear system with Jacobian 2292 $ endloop 2293 2294 Notes: 2295 Matrix must already be assemblied before calling this routine 2296 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2297 calling this routine. 2298 2299 .seealso: MatRetrieveValues() 2300 2301 @*/ 2302 int MatStoreValues(Mat mat) 2303 { 2304 int ierr,(*f)(Mat); 2305 2306 PetscFunctionBegin; 2307 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2308 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2309 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2310 2311 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2312 if (f) { 2313 ierr = (*f)(mat);CHKERRQ(ierr); 2314 } else { 2315 SETERRQ(1,"Wrong type of matrix to store values"); 2316 } 2317 PetscFunctionReturn(0); 2318 } 2319 2320 EXTERN_C_BEGIN 2321 #undef __FUNCT__ 2322 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2323 int MatRetrieveValues_SeqAIJ(Mat mat) 2324 { 2325 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2326 int nz = aij->i[mat->m],ierr; 2327 2328 PetscFunctionBegin; 2329 if (aij->nonew != 1) { 2330 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2331 } 2332 if (!aij->saved_values) { 2333 SETERRQ(1,"Must call MatStoreValues(A);first"); 2334 } 2335 2336 /* copy values over */ 2337 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2338 PetscFunctionReturn(0); 2339 } 2340 EXTERN_C_END 2341 2342 #undef __FUNCT__ 2343 #define __FUNCT__ "MatRetrieveValues" 2344 /*@ 2345 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2346 example, reuse of the linear part of a Jacobian, while recomputing the 2347 nonlinear portion. 2348 2349 Collect on Mat 2350 2351 Input Parameters: 2352 . mat - the matrix (currently on AIJ matrices support this option) 2353 2354 Level: advanced 2355 2356 .seealso: MatStoreValues() 2357 2358 @*/ 2359 int MatRetrieveValues(Mat mat) 2360 { 2361 int ierr,(*f)(Mat); 2362 2363 PetscFunctionBegin; 2364 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2365 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2366 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2367 2368 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2369 if (f) { 2370 ierr = (*f)(mat);CHKERRQ(ierr); 2371 } else { 2372 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2373 } 2374 PetscFunctionReturn(0); 2375 } 2376 2377 /* 2378 This allows SeqAIJ matrices to be passed to the matlab engine 2379 */ 2380 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2381 #include "engine.h" /* Matlab include file */ 2382 #include "mex.h" /* Matlab include file */ 2383 EXTERN_C_BEGIN 2384 #undef __FUNCT__ 2385 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2386 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 2387 { 2388 int ierr,i,*ai,*aj; 2389 Mat B = (Mat)obj; 2390 mxArray *mat; 2391 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2392 2393 PetscFunctionBegin; 2394 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2395 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2396 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2397 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2398 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2399 2400 /* Matlab indices start at 0 for sparse (what a surprise) */ 2401 2402 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2403 mxSetName(mat,obj->name); 2404 engPutArray((Engine *)mengine,mat); 2405 PetscFunctionReturn(0); 2406 } 2407 EXTERN_C_END 2408 2409 EXTERN_C_BEGIN 2410 #undef __FUNCT__ 2411 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2412 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 2413 { 2414 int ierr,ii; 2415 Mat mat = (Mat)obj; 2416 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2417 mxArray *mmat; 2418 2419 PetscFunctionBegin; 2420 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2421 2422 mmat = engGetArray((Engine *)mengine,obj->name); 2423 2424 aij->nz = (mxGetJc(mmat))[mat->m]; 2425 ierr = PetscMalloc(((size_t) aij->nz)*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2426 aij->j = (int*)(aij->a + aij->nz); 2427 aij->i = aij->j + aij->nz; 2428 aij->singlemalloc = PETSC_TRUE; 2429 aij->freedata = PETSC_TRUE; 2430 2431 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2432 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2433 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2434 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2435 2436 for (ii=0; ii<mat->m; ii++) { 2437 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2438 } 2439 2440 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2441 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2442 2443 PetscFunctionReturn(0); 2444 } 2445 EXTERN_C_END 2446 #endif 2447 2448 /* --------------------------------------------------------------------------------*/ 2449 #undef __FUNCT__ 2450 #define __FUNCT__ "MatCreateSeqAIJ" 2451 /*@C 2452 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2453 (the default parallel PETSc format). For good matrix assembly performance 2454 the user should preallocate the matrix storage by setting the parameter nz 2455 (or the array nnz). By setting these parameters accurately, performance 2456 during matrix assembly can be increased by more than a factor of 50. 2457 2458 Collective on MPI_Comm 2459 2460 Input Parameters: 2461 + comm - MPI communicator, set to PETSC_COMM_SELF 2462 . m - number of rows 2463 . n - number of columns 2464 . nz - number of nonzeros per row (same for all rows) 2465 - nnz - array containing the number of nonzeros in the various rows 2466 (possibly different for each row) or PETSC_NULL 2467 2468 Output Parameter: 2469 . A - the matrix 2470 2471 Notes: 2472 The AIJ format (also called the Yale sparse matrix format or 2473 compressed row storage), is fully compatible with standard Fortran 77 2474 storage. That is, the stored row and column indices can begin at 2475 either one (as in Fortran) or zero. See the users' manual for details. 2476 2477 Specify the preallocated storage with either nz or nnz (not both). 2478 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2479 allocation. For large problems you MUST preallocate memory or you 2480 will get TERRIBLE performance, see the users' manual chapter on matrices. 2481 2482 By default, this format uses inodes (identical nodes) when possible, to 2483 improve numerical efficiency of matrix-vector products and solves. We 2484 search for consecutive rows with the same nonzero structure, thereby 2485 reusing matrix information to achieve increased efficiency. 2486 2487 Options Database Keys: 2488 + -mat_aij_no_inode - Do not use inodes 2489 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2490 - -mat_aij_oneindex - Internally use indexing starting at 1 2491 rather than 0. Note that when calling MatSetValues(), 2492 the user still MUST index entries starting at 0! 2493 2494 Level: intermediate 2495 2496 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2497 2498 @*/ 2499 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,const int nnz[],Mat *A) 2500 { 2501 int ierr; 2502 2503 PetscFunctionBegin; 2504 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2505 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2506 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2507 PetscFunctionReturn(0); 2508 } 2509 2510 #define SKIP_ALLOCATION -4 2511 2512 #undef __FUNCT__ 2513 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2514 /*@C 2515 MatSeqAIJSetPreallocation - For good matrix assembly performance 2516 the user should preallocate the matrix storage by setting the parameter nz 2517 (or the array nnz). By setting these parameters accurately, performance 2518 during matrix assembly can be increased by more than a factor of 50. 2519 2520 Collective on MPI_Comm 2521 2522 Input Parameters: 2523 + comm - MPI communicator, set to PETSC_COMM_SELF 2524 . m - number of rows 2525 . n - number of columns 2526 . nz - number of nonzeros per row (same for all rows) 2527 - nnz - array containing the number of nonzeros in the various rows 2528 (possibly different for each row) or PETSC_NULL 2529 2530 Output Parameter: 2531 . A - the matrix 2532 2533 Notes: 2534 The AIJ format (also called the Yale sparse matrix format or 2535 compressed row storage), is fully compatible with standard Fortran 77 2536 storage. That is, the stored row and column indices can begin at 2537 either one (as in Fortran) or zero. See the users' manual for details. 2538 2539 Specify the preallocated storage with either nz or nnz (not both). 2540 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2541 allocation. For large problems you MUST preallocate memory or you 2542 will get TERRIBLE performance, see the users' manual chapter on matrices. 2543 2544 By default, this format uses inodes (identical nodes) when possible, to 2545 improve numerical efficiency of matrix-vector products and solves. We 2546 search for consecutive rows with the same nonzero structure, thereby 2547 reusing matrix information to achieve increased efficiency. 2548 2549 Options Database Keys: 2550 + -mat_aij_no_inode - Do not use inodes 2551 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2552 - -mat_aij_oneindex - Internally use indexing starting at 1 2553 rather than 0. Note that when calling MatSetValues(), 2554 the user still MUST index entries starting at 0! 2555 2556 Level: intermediate 2557 2558 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2559 2560 @*/ 2561 int MatSeqAIJSetPreallocation(Mat B,int nz,const int nnz[]) 2562 { 2563 int ierr,(*f)(Mat,int,const int[]); 2564 2565 PetscFunctionBegin; 2566 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2567 if (f) { 2568 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2569 } 2570 PetscFunctionReturn(0); 2571 } 2572 2573 EXTERN_C_BEGIN 2574 #undef __FUNCT__ 2575 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2576 int MatSeqAIJSetPreallocation_SeqAIJ(Mat B,int nz,int *nnz) 2577 { 2578 Mat_SeqAIJ *b; 2579 size_t len = 0; 2580 PetscTruth skipallocation = PETSC_FALSE; 2581 int i,ierr; 2582 2583 PetscFunctionBegin; 2584 2585 if (nz == SKIP_ALLOCATION) { 2586 skipallocation = PETSC_TRUE; 2587 nz = 0; 2588 } 2589 2590 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2591 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2592 if (nnz) { 2593 for (i=0; i<B->m; i++) { 2594 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2595 if (nnz[i] > B->n) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->n); 2596 } 2597 } 2598 2599 B->preallocated = PETSC_TRUE; 2600 b = (Mat_SeqAIJ*)B->data; 2601 2602 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2603 if (!nnz) { 2604 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2605 else if (nz <= 0) nz = 1; 2606 for (i=0; i<B->m; i++) b->imax[i] = nz; 2607 nz = nz*B->m; 2608 } else { 2609 nz = 0; 2610 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2611 } 2612 2613 if (!skipallocation) { 2614 /* allocate the matrix space */ 2615 len = ((size_t) nz)*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2616 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2617 b->j = (int*)(b->a + nz); 2618 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2619 b->i = b->j + nz; 2620 b->i[0] = 0; 2621 for (i=1; i<B->m+1; i++) { 2622 b->i[i] = b->i[i-1] + b->imax[i-1]; 2623 } 2624 b->singlemalloc = PETSC_TRUE; 2625 b->freedata = PETSC_TRUE; 2626 } else { 2627 b->freedata = PETSC_FALSE; 2628 } 2629 2630 /* b->ilen will count nonzeros in each row so far. */ 2631 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2632 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2633 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2634 2635 b->nz = 0; 2636 b->maxnz = nz; 2637 B->info.nz_unneeded = (double)b->maxnz; 2638 PetscFunctionReturn(0); 2639 } 2640 EXTERN_C_END 2641 2642 EXTERN int RegisterApplyPtAPRoutines_Private(Mat); 2643 2644 EXTERN_C_BEGIN 2645 EXTERN int MatConvert_SeqAIJ_SeqSBAIJ(Mat,MatType,Mat*); 2646 EXTERN int MatConvert_SeqAIJ_SeqBAIJ(Mat,MatType,Mat*); 2647 EXTERN int MatReorderForNonzeroDiagonal_SeqAIJ(Mat,PetscReal,IS,IS); 2648 EXTERN int MatAdjustForInodes_SeqAIJ(Mat,IS*,IS*); 2649 EXTERN int MatSeqAIJGetInodeSizes_SeqAIJ(Mat,int*,int*[],int*); 2650 EXTERN_C_END 2651 2652 EXTERN_C_BEGIN 2653 #undef __FUNCT__ 2654 #define __FUNCT__ "MatCreate_SeqAIJ" 2655 int MatCreate_SeqAIJ(Mat B) 2656 { 2657 Mat_SeqAIJ *b; 2658 int ierr,size; 2659 PetscTruth flg; 2660 2661 PetscFunctionBegin; 2662 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2663 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2664 2665 B->m = B->M = PetscMax(B->m,B->M); 2666 B->n = B->N = PetscMax(B->n,B->N); 2667 2668 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2669 B->data = (void*)b; 2670 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2671 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2672 B->factor = 0; 2673 B->lupivotthreshold = 1.0; 2674 B->mapping = 0; 2675 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2676 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2677 b->row = 0; 2678 b->col = 0; 2679 b->icol = 0; 2680 b->reallocs = 0; 2681 2682 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2683 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2684 2685 b->sorted = PETSC_FALSE; 2686 b->ignorezeroentries = PETSC_FALSE; 2687 b->roworiented = PETSC_TRUE; 2688 b->nonew = 0; 2689 b->diag = 0; 2690 b->solve_work = 0; 2691 B->spptr = 0; 2692 b->inode.use = PETSC_TRUE; 2693 b->inode.node_count = 0; 2694 b->inode.size = 0; 2695 b->inode.limit = 5; 2696 b->inode.max_limit = 5; 2697 b->saved_values = 0; 2698 b->idiag = 0; 2699 b->ssor = 0; 2700 b->keepzeroedrows = PETSC_FALSE; 2701 b->xtoy = 0; 2702 b->XtoY = 0; 2703 2704 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2705 2706 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2707 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2708 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2709 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2710 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2711 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2712 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2713 if (flg){ 2714 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2715 } 2716 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2717 "MatSeqAIJSetColumnIndices_SeqAIJ", 2718 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2719 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2720 "MatStoreValues_SeqAIJ", 2721 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2722 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2723 "MatRetrieveValues_SeqAIJ", 2724 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2725 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2726 "MatConvert_SeqAIJ_SeqSBAIJ", 2727 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2728 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2729 "MatConvert_SeqAIJ_SeqBAIJ", 2730 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2731 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsSymmetric_C", 2732 "MatIsSymmetric_SeqAIJ", 2733 MatIsSymmetric_SeqAIJ);CHKERRQ(ierr); 2734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2735 "MatSeqAIJSetPreallocation_SeqAIJ", 2736 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2737 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2738 "MatReorderForNonzeroDiagonal_SeqAIJ", 2739 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2740 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2741 "MatAdjustForInodes_SeqAIJ", 2742 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2743 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2744 "MatSeqAIJGetInodeSizes_SeqAIJ", 2745 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2746 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2747 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2748 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2749 #endif 2750 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2751 PetscFunctionReturn(0); 2752 } 2753 EXTERN_C_END 2754 2755 #undef __FUNCT__ 2756 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2757 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2758 { 2759 Mat C; 2760 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2761 int i,m = A->m,ierr; 2762 size_t len; 2763 2764 PetscFunctionBegin; 2765 *B = 0; 2766 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2767 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2768 c = (Mat_SeqAIJ*)C->data; 2769 2770 C->factor = A->factor; 2771 c->row = 0; 2772 c->col = 0; 2773 c->icol = 0; 2774 c->keepzeroedrows = a->keepzeroedrows; 2775 C->assembled = PETSC_TRUE; 2776 2777 C->M = A->m; 2778 C->N = A->n; 2779 2780 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2781 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2782 for (i=0; i<m; i++) { 2783 c->imax[i] = a->imax[i]; 2784 c->ilen[i] = a->ilen[i]; 2785 } 2786 2787 /* allocate the matrix space */ 2788 c->singlemalloc = PETSC_TRUE; 2789 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2790 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2791 c->j = (int*)(c->a + a->i[m] ); 2792 c->i = c->j + a->i[m]; 2793 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2794 if (m > 0) { 2795 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(int));CHKERRQ(ierr); 2796 if (cpvalues == MAT_COPY_VALUES) { 2797 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2798 } else { 2799 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2800 } 2801 } 2802 2803 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2804 c->sorted = a->sorted; 2805 c->roworiented = a->roworiented; 2806 c->nonew = a->nonew; 2807 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2808 c->saved_values = 0; 2809 c->idiag = 0; 2810 c->ssor = 0; 2811 c->ignorezeroentries = a->ignorezeroentries; 2812 c->freedata = PETSC_TRUE; 2813 2814 if (a->diag) { 2815 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2816 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2817 for (i=0; i<m; i++) { 2818 c->diag[i] = a->diag[i]; 2819 } 2820 } else c->diag = 0; 2821 c->inode.use = a->inode.use; 2822 c->inode.limit = a->inode.limit; 2823 c->inode.max_limit = a->inode.max_limit; 2824 if (a->inode.size){ 2825 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2826 c->inode.node_count = a->inode.node_count; 2827 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2828 } else { 2829 c->inode.size = 0; 2830 c->inode.node_count = 0; 2831 } 2832 c->nz = a->nz; 2833 c->maxnz = a->maxnz; 2834 c->solve_work = 0; 2835 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2836 C->preallocated = PETSC_TRUE; 2837 2838 *B = C; 2839 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2840 PetscFunctionReturn(0); 2841 } 2842 2843 EXTERN_C_BEGIN 2844 #undef __FUNCT__ 2845 #define __FUNCT__ "MatLoad_SeqAIJ" 2846 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2847 { 2848 Mat_SeqAIJ *a; 2849 Mat B; 2850 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N; 2851 MPI_Comm comm; 2852 #if defined(PETSC_HAVE_UMFPACK) || defined(PETSC_HAVE_MUMPS) 2853 PetscTruth flag; 2854 #endif 2855 2856 PetscFunctionBegin; 2857 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2858 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2859 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2860 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2861 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2862 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2863 M = header[1]; N = header[2]; nz = header[3]; 2864 2865 if (nz < 0) { 2866 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2867 } 2868 2869 /* read in row lengths */ 2870 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2871 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2872 2873 /* create our matrix */ 2874 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2875 ierr = MatSetType(B,type);CHKERRQ(ierr); 2876 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2877 a = (Mat_SeqAIJ*)B->data; 2878 2879 /* read in column indices and adjust for Fortran indexing*/ 2880 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2881 2882 /* read in nonzero values */ 2883 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2884 2885 /* set matrix "i" values */ 2886 a->i[0] = 0; 2887 for (i=1; i<= M; i++) { 2888 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2889 a->ilen[i-1] = rowlengths[i-1]; 2890 } 2891 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2892 2893 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2894 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2895 #if defined(PETSC_HAVE_UMFPACK) 2896 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_umfpack",&flag);CHKERRQ(ierr); 2897 if (flag) { ierr = MatUseUMFPACK_SeqAIJ(B);CHKERRQ(ierr); } 2898 #endif 2899 #if defined(PETSC_HAVE_MUMPS) 2900 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_mumps",&flag);CHKERRQ(ierr); 2901 if (flag) { ierr = MatUseMUMPS_MPIAIJ(B);CHKERRQ(ierr); } 2902 #endif 2903 *A = B; 2904 PetscFunctionReturn(0); 2905 } 2906 EXTERN_C_END 2907 2908 #undef __FUNCT__ 2909 #define __FUNCT__ "MatEqual_SeqAIJ" 2910 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2911 { 2912 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2913 int ierr; 2914 2915 PetscFunctionBegin; 2916 /* If the matrix dimensions are not equal,or no of nonzeros */ 2917 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2918 *flg = PETSC_FALSE; 2919 PetscFunctionReturn(0); 2920 } 2921 2922 /* if the a->i are the same */ 2923 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2924 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2925 2926 /* if a->j are the same */ 2927 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2928 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2929 2930 /* if a->a are the same */ 2931 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2932 2933 PetscFunctionReturn(0); 2934 2935 } 2936 2937 #undef __FUNCT__ 2938 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2939 /*@C 2940 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2941 provided by the user. 2942 2943 Coolective on MPI_Comm 2944 2945 Input Parameters: 2946 + comm - must be an MPI communicator of size 1 2947 . m - number of rows 2948 . n - number of columns 2949 . i - row indices 2950 . j - column indices 2951 - a - matrix values 2952 2953 Output Parameter: 2954 . mat - the matrix 2955 2956 Level: intermediate 2957 2958 Notes: 2959 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2960 once the matrix is destroyed 2961 2962 You cannot set new nonzero locations into this matrix, that will generate an error. 2963 2964 The i and j indices are 0 based 2965 2966 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2967 2968 @*/ 2969 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2970 { 2971 int ierr,ii; 2972 Mat_SeqAIJ *aij; 2973 2974 PetscFunctionBegin; 2975 ierr = MatCreateSeqAIJ(comm,m,n,SKIP_ALLOCATION,0,mat);CHKERRQ(ierr); 2976 aij = (Mat_SeqAIJ*)(*mat)->data; 2977 2978 if (i[0] != 0) { 2979 SETERRQ(1,"i (row indices) must start with 0"); 2980 } 2981 aij->i = i; 2982 aij->j = j; 2983 aij->a = a; 2984 aij->singlemalloc = PETSC_FALSE; 2985 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2986 aij->freedata = PETSC_FALSE; 2987 2988 for (ii=0; ii<m; ii++) { 2989 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2990 #if defined(PETSC_USE_BOPT_g) 2991 if (i[ii+1] - i[ii] < 0) SETERRQ2(1,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2992 #endif 2993 } 2994 #if defined(PETSC_USE_BOPT_g) 2995 for (ii=0; ii<aij->i[m]; ii++) { 2996 if (j[ii] < 0) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2997 if (j[ii] > n - 1) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2998 } 2999 #endif 3000 3001 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3002 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3003 PetscFunctionReturn(0); 3004 } 3005 3006 #undef __FUNCT__ 3007 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3008 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3009 { 3010 int ierr; 3011 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3012 3013 PetscFunctionBegin; 3014 if (coloring->ctype == IS_COLORING_LOCAL) { 3015 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3016 a->coloring = coloring; 3017 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3018 int i,*larray; 3019 ISColoring ocoloring; 3020 ISColoringValue *colors; 3021 3022 /* set coloring for diagonal portion */ 3023 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 3024 for (i=0; i<A->n; i++) { 3025 larray[i] = i; 3026 } 3027 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3028 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3029 for (i=0; i<A->n; i++) { 3030 colors[i] = coloring->colors[larray[i]]; 3031 } 3032 ierr = PetscFree(larray);CHKERRQ(ierr); 3033 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3034 a->coloring = ocoloring; 3035 } 3036 PetscFunctionReturn(0); 3037 } 3038 3039 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 3040 EXTERN_C_BEGIN 3041 #include "adic/ad_utils.h" 3042 EXTERN_C_END 3043 3044 #undef __FUNCT__ 3045 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3046 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3047 { 3048 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3049 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3050 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3051 ISColoringValue *color; 3052 3053 PetscFunctionBegin; 3054 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3055 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3056 color = a->coloring->colors; 3057 /* loop over rows */ 3058 for (i=0; i<m; i++) { 3059 nz = ii[i+1] - ii[i]; 3060 /* loop over columns putting computed value into matrix */ 3061 for (j=0; j<nz; j++) { 3062 *v++ = values[color[*jj++]]; 3063 } 3064 values += nlen; /* jump to next row of derivatives */ 3065 } 3066 PetscFunctionReturn(0); 3067 } 3068 3069 #else 3070 3071 #undef __FUNCT__ 3072 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3073 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3074 { 3075 PetscFunctionBegin; 3076 SETERRQ(1,"PETSc installed without ADIC"); 3077 } 3078 3079 #endif 3080 3081 #undef __FUNCT__ 3082 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3083 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 3084 { 3085 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3086 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3087 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3088 ISColoringValue *color; 3089 3090 PetscFunctionBegin; 3091 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3092 color = a->coloring->colors; 3093 /* loop over rows */ 3094 for (i=0; i<m; i++) { 3095 nz = ii[i+1] - ii[i]; 3096 /* loop over columns putting computed value into matrix */ 3097 for (j=0; j<nz; j++) { 3098 *v++ = values[color[*jj++]]; 3099 } 3100 values += nl; /* jump to next row of derivatives */ 3101 } 3102 PetscFunctionReturn(0); 3103 } 3104 3105