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