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