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 = VecGetArray(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 = VecRestoreArray(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 = VecGetArray(xx,&x);CHKERRQ(ierr); 826 ierr = VecGetArray(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 = VecRestoreArray(xx,&x);CHKERRQ(ierr); 841 ierr = VecRestoreArray(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 = VecGetArray(xx,&x);CHKERRQ(ierr); 877 ierr = VecGetArray(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 = VecRestoreArray(xx,&x);CHKERRQ(ierr); 896 ierr = VecRestoreArray(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 = VecGetArray(xx,&x);CHKERRQ(ierr); 914 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 915 if (zz != yy) { 916 ierr = VecGetArray(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 = VecRestoreArray(xx,&x);CHKERRQ(ierr); 939 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 940 if (zz != yy) { 941 ierr = VecRestoreArray(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 = VecGetArray(xx,&x);CHKERRQ(ierr); 1038 if (xx != bb) { 1039 ierr = VecGetArray(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 = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1059 if (bb != xx) {ierr = VecRestoreArray(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 = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1119 if (bb != xx) {ierr = VecRestoreArray(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 = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1198 if (bb != xx) {ierr = VecRestoreArray(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 = MatCreateSeqAIJ(A->comm,A->n,m,0,col,&C);CHKERRQ(ierr); 1377 ierr = PetscFree(col);CHKERRQ(ierr); 1378 for (i=0; i<m; i++) { 1379 len = ai[i+1]-ai[i]; 1380 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1381 array += len; 1382 aj += len; 1383 } 1384 1385 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1386 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1387 1388 if (B) { 1389 *B = C; 1390 } else { 1391 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1392 } 1393 PetscFunctionReturn(0); 1394 } 1395 1396 EXTERN_C_BEGIN 1397 #undef __FUNCT__ 1398 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1399 int MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscTruth *f) 1400 { 1401 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1402 int *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1403 int ma,na,mb,nb, i,ierr; 1404 1405 PetscFunctionBegin; 1406 bij = (Mat_SeqAIJ *) B->data; 1407 1408 ierr = MatGetSize(A,&ma,&na); CHKERRQ(ierr); 1409 ierr = MatGetSize(B,&mb,&nb); CHKERRQ(ierr); 1410 if (ma!=nb || na!=mb) 1411 SETERRQ(1,"Incompatible A/B sizes for symmetry test"); 1412 aii = aij->i; bii = bij->i; 1413 adx = aij->j; bdx = bij->j; 1414 va = aij->a; vb = bij->a; 1415 ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr); 1416 ierr = PetscMalloc(mb*sizeof(int),&bptr); CHKERRQ(ierr); 1417 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1418 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1419 1420 *f = PETSC_TRUE; 1421 for (i=0; i<ma; i++) { 1422 /*printf("row %d spans %d--%d; we start @ %d\n", 1423 i,idx[ii[i]],idx[ii[i+1]-1],idx[aptr[i]]);*/ 1424 while (aptr[i]<aii[i+1]) { 1425 int idc,idr; PetscScalar vc,vr; 1426 /* column/row index/value */ 1427 idc = adx[aptr[i]]; idr = bdx[bptr[idc]]; 1428 vc = va[aptr[i]]; vr = vb[bptr[idc]]; 1429 /*printf("comparing %d: (%d,%d)=%e to (%d,%d)=%e\n", 1430 aptr[i],i,idc,vc,idc,idr,vr);*/ 1431 if (i!=idr || vc!=vr) { 1432 *f = PETSC_FALSE; goto done; 1433 } else { 1434 aptr[i]++; if (B || i!=idc) bptr[idc]++; 1435 } 1436 } 1437 } 1438 done: 1439 ierr = PetscFree(aptr); CHKERRQ(ierr); 1440 if (B) { 1441 ierr = PetscFree(bptr); CHKERRQ(ierr); 1442 } 1443 1444 PetscFunctionReturn(0); 1445 } 1446 EXTERN_C_END 1447 1448 #undef __FUNCT__ 1449 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1450 int MatIsSymmetric_SeqAIJ(Mat A,PetscTruth *f) 1451 { 1452 int ierr; 1453 PetscFunctionBegin; 1454 ierr = MatIsTranspose_SeqAIJ(A,A,f); CHKERRQ(ierr); 1455 PetscFunctionReturn(0); 1456 } 1457 1458 #undef __FUNCT__ 1459 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1460 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1461 { 1462 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1463 PetscScalar *l,*r,x,*v; 1464 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1465 1466 PetscFunctionBegin; 1467 if (ll) { 1468 /* The local size is used so that VecMPI can be passed to this routine 1469 by MatDiagonalScale_MPIAIJ */ 1470 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1471 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1472 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1473 v = a->a; 1474 for (i=0; i<m; i++) { 1475 x = l[i]; 1476 M = a->i[i+1] - a->i[i]; 1477 for (j=0; j<M; j++) { (*v++) *= x;} 1478 } 1479 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1480 PetscLogFlops(nz); 1481 } 1482 if (rr) { 1483 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1484 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1485 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1486 v = a->a; jj = a->j; 1487 for (i=0; i<nz; i++) { 1488 (*v++) *= r[*jj++]; 1489 } 1490 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1491 PetscLogFlops(nz); 1492 } 1493 PetscFunctionReturn(0); 1494 } 1495 1496 #undef __FUNCT__ 1497 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1498 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1499 { 1500 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1501 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1502 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1503 int *irow,*icol,nrows,ncols; 1504 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1505 PetscScalar *a_new,*mat_a; 1506 Mat C; 1507 PetscTruth stride; 1508 1509 PetscFunctionBegin; 1510 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1511 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1512 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1513 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1514 1515 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1516 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1517 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1518 1519 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1520 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1521 if (stride && step == 1) { 1522 /* special case of contiguous rows */ 1523 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1524 starts = lens + nrows; 1525 /* loop over new rows determining lens and starting points */ 1526 for (i=0; i<nrows; i++) { 1527 kstart = ai[irow[i]]; 1528 kend = kstart + ailen[irow[i]]; 1529 for (k=kstart; k<kend; k++) { 1530 if (aj[k] >= first) { 1531 starts[i] = k; 1532 break; 1533 } 1534 } 1535 sum = 0; 1536 while (k < kend) { 1537 if (aj[k++] >= first+ncols) break; 1538 sum++; 1539 } 1540 lens[i] = sum; 1541 } 1542 /* create submatrix */ 1543 if (scall == MAT_REUSE_MATRIX) { 1544 int n_cols,n_rows; 1545 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1546 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1547 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1548 C = *B; 1549 } else { 1550 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1551 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1552 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1553 } 1554 c = (Mat_SeqAIJ*)C->data; 1555 1556 /* loop over rows inserting into submatrix */ 1557 a_new = c->a; 1558 j_new = c->j; 1559 i_new = c->i; 1560 1561 for (i=0; i<nrows; i++) { 1562 ii = starts[i]; 1563 lensi = lens[i]; 1564 for (k=0; k<lensi; k++) { 1565 *j_new++ = aj[ii+k] - first; 1566 } 1567 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1568 a_new += lensi; 1569 i_new[i+1] = i_new[i] + lensi; 1570 c->ilen[i] = lensi; 1571 } 1572 ierr = PetscFree(lens);CHKERRQ(ierr); 1573 } else { 1574 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1575 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1576 1577 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1578 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1579 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1580 /* determine lens of each row */ 1581 for (i=0; i<nrows; i++) { 1582 kstart = ai[irow[i]]; 1583 kend = kstart + a->ilen[irow[i]]; 1584 lens[i] = 0; 1585 for (k=kstart; k<kend; k++) { 1586 if (smap[aj[k]]) { 1587 lens[i]++; 1588 } 1589 } 1590 } 1591 /* Create and fill new matrix */ 1592 if (scall == MAT_REUSE_MATRIX) { 1593 PetscTruth equal; 1594 1595 c = (Mat_SeqAIJ *)((*B)->data); 1596 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1597 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1598 if (!equal) { 1599 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1600 } 1601 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1602 C = *B; 1603 } else { 1604 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1605 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1606 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1607 } 1608 c = (Mat_SeqAIJ *)(C->data); 1609 for (i=0; i<nrows; i++) { 1610 row = irow[i]; 1611 kstart = ai[row]; 1612 kend = kstart + a->ilen[row]; 1613 mat_i = c->i[i]; 1614 mat_j = c->j + mat_i; 1615 mat_a = c->a + mat_i; 1616 mat_ilen = c->ilen + i; 1617 for (k=kstart; k<kend; k++) { 1618 if ((tcol=smap[a->j[k]])) { 1619 *mat_j++ = tcol - 1; 1620 *mat_a++ = a->a[k]; 1621 (*mat_ilen)++; 1622 1623 } 1624 } 1625 } 1626 /* Free work space */ 1627 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1628 ierr = PetscFree(smap);CHKERRQ(ierr); 1629 ierr = PetscFree(lens);CHKERRQ(ierr); 1630 } 1631 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1632 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1633 1634 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1635 *B = C; 1636 PetscFunctionReturn(0); 1637 } 1638 1639 /* 1640 */ 1641 #undef __FUNCT__ 1642 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1643 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1644 { 1645 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1646 int ierr; 1647 Mat outA; 1648 PetscTruth row_identity,col_identity; 1649 1650 PetscFunctionBegin; 1651 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1652 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1653 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1654 if (!row_identity || !col_identity) { 1655 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1656 } 1657 1658 outA = inA; 1659 inA->factor = FACTOR_LU; 1660 a->row = row; 1661 a->col = col; 1662 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1663 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1664 1665 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1666 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1667 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1668 PetscLogObjectParent(inA,a->icol); 1669 1670 if (!a->solve_work) { /* this matrix may have been factored before */ 1671 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1672 } 1673 1674 if (!a->diag) { 1675 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1676 } 1677 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1678 PetscFunctionReturn(0); 1679 } 1680 1681 #include "petscblaslapack.h" 1682 #undef __FUNCT__ 1683 #define __FUNCT__ "MatScale_SeqAIJ" 1684 int MatScale_SeqAIJ(const PetscScalar *alpha,Mat inA) 1685 { 1686 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1687 int one = 1; 1688 1689 PetscFunctionBegin; 1690 BLscal_(&a->nz,(PetscScalar*)alpha,a->a,&one); 1691 PetscLogFlops(a->nz); 1692 PetscFunctionReturn(0); 1693 } 1694 1695 #undef __FUNCT__ 1696 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1697 int MatGetSubMatrices_SeqAIJ(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1698 { 1699 int ierr,i; 1700 1701 PetscFunctionBegin; 1702 if (scall == MAT_INITIAL_MATRIX) { 1703 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1704 } 1705 1706 for (i=0; i<n; i++) { 1707 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1708 } 1709 PetscFunctionReturn(0); 1710 } 1711 1712 #undef __FUNCT__ 1713 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1714 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1715 { 1716 PetscFunctionBegin; 1717 *bs = 1; 1718 PetscFunctionReturn(0); 1719 } 1720 1721 #undef __FUNCT__ 1722 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1723 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS is[],int ov) 1724 { 1725 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1726 int row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1727 int start,end,*ai,*aj; 1728 PetscBT table; 1729 1730 PetscFunctionBegin; 1731 m = A->m; 1732 ai = a->i; 1733 aj = a->j; 1734 1735 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1736 1737 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1738 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1739 1740 for (i=0; i<is_max; i++) { 1741 /* Initialize the two local arrays */ 1742 isz = 0; 1743 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1744 1745 /* Extract the indices, assume there can be duplicate entries */ 1746 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1747 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1748 1749 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1750 for (j=0; j<n ; ++j){ 1751 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1752 } 1753 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1754 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1755 1756 k = 0; 1757 for (j=0; j<ov; j++){ /* for each overlap */ 1758 n = isz; 1759 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1760 row = nidx[k]; 1761 start = ai[row]; 1762 end = ai[row+1]; 1763 for (l = start; l<end ; l++){ 1764 val = aj[l] ; 1765 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1766 } 1767 } 1768 } 1769 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1770 } 1771 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1772 ierr = PetscFree(nidx);CHKERRQ(ierr); 1773 PetscFunctionReturn(0); 1774 } 1775 1776 /* -------------------------------------------------------------- */ 1777 #undef __FUNCT__ 1778 #define __FUNCT__ "MatPermute_SeqAIJ" 1779 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1780 { 1781 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1782 PetscScalar *vwork; 1783 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1784 int *row,*col,*cnew,j,*lens; 1785 IS icolp,irowp; 1786 1787 PetscFunctionBegin; 1788 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1789 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1790 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1791 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1792 1793 /* determine lengths of permuted rows */ 1794 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1795 for (i=0; i<m; i++) { 1796 lens[row[i]] = a->i[i+1] - a->i[i]; 1797 } 1798 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1799 ierr = PetscFree(lens);CHKERRQ(ierr); 1800 1801 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1802 for (i=0; i<m; i++) { 1803 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1804 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1805 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1806 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1807 } 1808 ierr = PetscFree(cnew);CHKERRQ(ierr); 1809 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1810 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1811 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1812 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1813 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1814 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1815 PetscFunctionReturn(0); 1816 } 1817 1818 #undef __FUNCT__ 1819 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1820 int MatPrintHelp_SeqAIJ(Mat A) 1821 { 1822 static PetscTruth called = PETSC_FALSE; 1823 MPI_Comm comm = A->comm; 1824 int ierr; 1825 1826 PetscFunctionBegin; 1827 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1828 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1829 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1830 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1831 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1832 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1833 #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1834 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1835 #endif 1836 PetscFunctionReturn(0); 1837 } 1838 1839 #undef __FUNCT__ 1840 #define __FUNCT__ "MatCopy_SeqAIJ" 1841 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1842 { 1843 int ierr; 1844 1845 PetscFunctionBegin; 1846 /* If the two matrices have the same copy implementation, use fast copy. */ 1847 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1848 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1849 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1850 1851 if (a->i[A->m] != b->i[B->m]) { 1852 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1853 } 1854 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1855 } else { 1856 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1857 } 1858 PetscFunctionReturn(0); 1859 } 1860 1861 #undef __FUNCT__ 1862 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1863 int MatSetUpPreallocation_SeqAIJ(Mat A) 1864 { 1865 int ierr; 1866 1867 PetscFunctionBegin; 1868 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1869 PetscFunctionReturn(0); 1870 } 1871 1872 #undef __FUNCT__ 1873 #define __FUNCT__ "MatGetArray_SeqAIJ" 1874 int MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1875 { 1876 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1877 PetscFunctionBegin; 1878 *array = a->a; 1879 PetscFunctionReturn(0); 1880 } 1881 1882 #undef __FUNCT__ 1883 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1884 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1885 { 1886 PetscFunctionBegin; 1887 PetscFunctionReturn(0); 1888 } 1889 1890 #undef __FUNCT__ 1891 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1892 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1893 { 1894 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1895 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1896 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1897 PetscScalar *vscale_array; 1898 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1899 Vec w1,w2,w3; 1900 void *fctx = coloring->fctx; 1901 PetscTruth flg; 1902 1903 PetscFunctionBegin; 1904 if (!coloring->w1) { 1905 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1906 PetscLogObjectParent(coloring,coloring->w1); 1907 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1908 PetscLogObjectParent(coloring,coloring->w2); 1909 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1910 PetscLogObjectParent(coloring,coloring->w3); 1911 } 1912 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1913 1914 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1915 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1916 if (flg) { 1917 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1918 } else { 1919 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1920 } 1921 1922 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1923 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1924 1925 /* 1926 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1927 coloring->F for the coarser grids from the finest 1928 */ 1929 if (coloring->F) { 1930 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1931 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1932 if (m1 != m2) { 1933 coloring->F = 0; 1934 } 1935 } 1936 1937 if (coloring->F) { 1938 w1 = coloring->F; 1939 coloring->F = 0; 1940 } else { 1941 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1942 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1943 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1944 } 1945 1946 /* 1947 Compute all the scale factors and share with other processors 1948 */ 1949 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1950 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1951 for (k=0; k<coloring->ncolors; k++) { 1952 /* 1953 Loop over each column associated with color adding the 1954 perturbation to the vector w3. 1955 */ 1956 for (l=0; l<coloring->ncolumns[k]; l++) { 1957 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1958 dx = xx[col]; 1959 if (dx == 0.0) dx = 1.0; 1960 #if !defined(PETSC_USE_COMPLEX) 1961 if (dx < umin && dx >= 0.0) dx = umin; 1962 else if (dx < 0.0 && dx > -umin) dx = -umin; 1963 #else 1964 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1965 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1966 #endif 1967 dx *= epsilon; 1968 vscale_array[col] = 1.0/dx; 1969 } 1970 } 1971 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1972 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1973 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1974 1975 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1976 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 1977 1978 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 1979 else vscaleforrow = coloring->columnsforrow; 1980 1981 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1982 /* 1983 Loop over each color 1984 */ 1985 for (k=0; k<coloring->ncolors; k++) { 1986 coloring->currentcolor = k; 1987 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 1988 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 1989 /* 1990 Loop over each column associated with color adding the 1991 perturbation to the vector w3. 1992 */ 1993 for (l=0; l<coloring->ncolumns[k]; l++) { 1994 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1995 dx = xx[col]; 1996 if (dx == 0.0) dx = 1.0; 1997 #if !defined(PETSC_USE_COMPLEX) 1998 if (dx < umin && dx >= 0.0) dx = umin; 1999 else if (dx < 0.0 && dx > -umin) dx = -umin; 2000 #else 2001 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2002 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2003 #endif 2004 dx *= epsilon; 2005 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2006 w3_array[col] += dx; 2007 } 2008 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2009 2010 /* 2011 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2012 */ 2013 2014 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2015 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2016 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2017 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2018 2019 /* 2020 Loop over rows of vector, putting results into Jacobian matrix 2021 */ 2022 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2023 for (l=0; l<coloring->nrows[k]; l++) { 2024 row = coloring->rows[k][l]; 2025 col = coloring->columnsforrow[k][l]; 2026 y[row] *= vscale_array[vscaleforrow[k][l]]; 2027 srow = row + start; 2028 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2029 } 2030 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2031 } 2032 coloring->currentcolor = k; 2033 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2034 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2035 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2036 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2037 PetscFunctionReturn(0); 2038 } 2039 2040 #include "petscblaslapack.h" 2041 #undef __FUNCT__ 2042 #define __FUNCT__ "MatAXPY_SeqAIJ" 2043 int MatAXPY_SeqAIJ(const PetscScalar a[],Mat X,Mat Y,MatStructure str) 2044 { 2045 int ierr,one=1,i; 2046 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2047 2048 PetscFunctionBegin; 2049 if (str == SAME_NONZERO_PATTERN) { 2050 BLaxpy_(&x->nz,(PetscScalar*)a,x->a,&one,y->a,&one); 2051 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2052 if (y->xtoy && y->XtoY != X) { 2053 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2054 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2055 } 2056 if (!y->xtoy) { /* get xtoy */ 2057 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2058 y->XtoY = X; 2059 } 2060 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2061 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2062 } else { 2063 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2064 } 2065 PetscFunctionReturn(0); 2066 } 2067 2068 /* -------------------------------------------------------------------*/ 2069 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2070 MatGetRow_SeqAIJ, 2071 MatRestoreRow_SeqAIJ, 2072 MatMult_SeqAIJ, 2073 /* 4*/ MatMultAdd_SeqAIJ, 2074 MatMultTranspose_SeqAIJ, 2075 MatMultTransposeAdd_SeqAIJ, 2076 MatSolve_SeqAIJ, 2077 MatSolveAdd_SeqAIJ, 2078 MatSolveTranspose_SeqAIJ, 2079 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2080 MatLUFactor_SeqAIJ, 2081 0, 2082 MatRelax_SeqAIJ, 2083 MatTranspose_SeqAIJ, 2084 /*15*/ MatGetInfo_SeqAIJ, 2085 MatEqual_SeqAIJ, 2086 MatGetDiagonal_SeqAIJ, 2087 MatDiagonalScale_SeqAIJ, 2088 MatNorm_SeqAIJ, 2089 /*20*/ 0, 2090 MatAssemblyEnd_SeqAIJ, 2091 MatCompress_SeqAIJ, 2092 MatSetOption_SeqAIJ, 2093 MatZeroEntries_SeqAIJ, 2094 /*25*/ MatZeroRows_SeqAIJ, 2095 MatLUFactorSymbolic_SeqAIJ, 2096 MatLUFactorNumeric_SeqAIJ, 2097 MatCholeskyFactorSymbolic_SeqAIJ, 2098 MatCholeskyFactorNumeric_SeqAIJ, 2099 /*30*/ MatSetUpPreallocation_SeqAIJ, 2100 MatILUFactorSymbolic_SeqAIJ, 2101 MatICCFactorSymbolic_SeqAIJ, 2102 MatGetArray_SeqAIJ, 2103 MatRestoreArray_SeqAIJ, 2104 /*35*/ MatDuplicate_SeqAIJ, 2105 0, 2106 0, 2107 MatILUFactor_SeqAIJ, 2108 0, 2109 /*40*/ MatAXPY_SeqAIJ, 2110 MatGetSubMatrices_SeqAIJ, 2111 MatIncreaseOverlap_SeqAIJ, 2112 MatGetValues_SeqAIJ, 2113 MatCopy_SeqAIJ, 2114 /*45*/ MatPrintHelp_SeqAIJ, 2115 MatScale_SeqAIJ, 2116 0, 2117 0, 2118 MatILUDTFactor_SeqAIJ, 2119 /*50*/ MatGetBlockSize_SeqAIJ, 2120 MatGetRowIJ_SeqAIJ, 2121 MatRestoreRowIJ_SeqAIJ, 2122 MatGetColumnIJ_SeqAIJ, 2123 MatRestoreColumnIJ_SeqAIJ, 2124 /*55*/ MatFDColoringCreate_SeqAIJ, 2125 0, 2126 0, 2127 MatPermute_SeqAIJ, 2128 0, 2129 /*60*/ 0, 2130 MatDestroy_SeqAIJ, 2131 MatView_SeqAIJ, 2132 MatGetPetscMaps_Petsc, 2133 0, 2134 /*65*/ 0, 2135 0, 2136 0, 2137 0, 2138 0, 2139 /*70*/ 0, 2140 0, 2141 MatSetColoring_SeqAIJ, 2142 MatSetValuesAdic_SeqAIJ, 2143 MatSetValuesAdifor_SeqAIJ, 2144 /*75*/ MatFDColoringApply_SeqAIJ, 2145 0, 2146 0, 2147 0, 2148 0, 2149 /*80*/ 0, 2150 0, 2151 0, 2152 0, 2153 /*85*/ MatLoad_SeqAIJ, 2154 MatIsSymmetric_SeqAIJ, 2155 }; 2156 2157 EXTERN_C_BEGIN 2158 #undef __FUNCT__ 2159 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2160 2161 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2162 { 2163 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2164 int i,nz,n; 2165 2166 PetscFunctionBegin; 2167 2168 nz = aij->maxnz; 2169 n = mat->n; 2170 for (i=0; i<nz; i++) { 2171 aij->j[i] = indices[i]; 2172 } 2173 aij->nz = nz; 2174 for (i=0; i<n; i++) { 2175 aij->ilen[i] = aij->imax[i]; 2176 } 2177 2178 PetscFunctionReturn(0); 2179 } 2180 EXTERN_C_END 2181 2182 #undef __FUNCT__ 2183 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2184 /*@ 2185 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2186 in the matrix. 2187 2188 Input Parameters: 2189 + mat - the SeqAIJ matrix 2190 - indices - the column indices 2191 2192 Level: advanced 2193 2194 Notes: 2195 This can be called if you have precomputed the nonzero structure of the 2196 matrix and want to provide it to the matrix object to improve the performance 2197 of the MatSetValues() operation. 2198 2199 You MUST have set the correct numbers of nonzeros per row in the call to 2200 MatCreateSeqAIJ(). 2201 2202 MUST be called before any calls to MatSetValues(); 2203 2204 The indices should start with zero, not one. 2205 2206 @*/ 2207 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2208 { 2209 int ierr,(*f)(Mat,int *); 2210 2211 PetscFunctionBegin; 2212 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2213 PetscValidPointer(indices,2); 2214 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2215 if (f) { 2216 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2217 } else { 2218 SETERRQ(1,"Wrong type of matrix to set column indices"); 2219 } 2220 PetscFunctionReturn(0); 2221 } 2222 2223 /* ----------------------------------------------------------------------------------------*/ 2224 2225 EXTERN_C_BEGIN 2226 #undef __FUNCT__ 2227 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2228 int MatStoreValues_SeqAIJ(Mat mat) 2229 { 2230 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2231 size_t nz = aij->i[mat->m],ierr; 2232 2233 PetscFunctionBegin; 2234 if (aij->nonew != 1) { 2235 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2236 } 2237 2238 /* allocate space for values if not already there */ 2239 if (!aij->saved_values) { 2240 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2241 } 2242 2243 /* copy values over */ 2244 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2245 PetscFunctionReturn(0); 2246 } 2247 EXTERN_C_END 2248 2249 #undef __FUNCT__ 2250 #define __FUNCT__ "MatStoreValues" 2251 /*@ 2252 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2253 example, reuse of the linear part of a Jacobian, while recomputing the 2254 nonlinear portion. 2255 2256 Collect on Mat 2257 2258 Input Parameters: 2259 . mat - the matrix (currently on AIJ matrices support this option) 2260 2261 Level: advanced 2262 2263 Common Usage, with SNESSolve(): 2264 $ Create Jacobian matrix 2265 $ Set linear terms into matrix 2266 $ Apply boundary conditions to matrix, at this time matrix must have 2267 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2268 $ boundary conditions again will not change the nonzero structure 2269 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2270 $ ierr = MatStoreValues(mat); 2271 $ Call SNESSetJacobian() with matrix 2272 $ In your Jacobian routine 2273 $ ierr = MatRetrieveValues(mat); 2274 $ Set nonlinear terms in matrix 2275 2276 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2277 $ // build linear portion of Jacobian 2278 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2279 $ ierr = MatStoreValues(mat); 2280 $ loop over nonlinear iterations 2281 $ ierr = MatRetrieveValues(mat); 2282 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2283 $ // call MatAssemblyBegin/End() on matrix 2284 $ Solve linear system with Jacobian 2285 $ endloop 2286 2287 Notes: 2288 Matrix must already be assemblied before calling this routine 2289 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2290 calling this routine. 2291 2292 .seealso: MatRetrieveValues() 2293 2294 @*/ 2295 int MatStoreValues(Mat mat) 2296 { 2297 int ierr,(*f)(Mat); 2298 2299 PetscFunctionBegin; 2300 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2301 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2302 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2303 2304 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2305 if (f) { 2306 ierr = (*f)(mat);CHKERRQ(ierr); 2307 } else { 2308 SETERRQ(1,"Wrong type of matrix to store values"); 2309 } 2310 PetscFunctionReturn(0); 2311 } 2312 2313 EXTERN_C_BEGIN 2314 #undef __FUNCT__ 2315 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2316 int MatRetrieveValues_SeqAIJ(Mat mat) 2317 { 2318 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2319 int nz = aij->i[mat->m],ierr; 2320 2321 PetscFunctionBegin; 2322 if (aij->nonew != 1) { 2323 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2324 } 2325 if (!aij->saved_values) { 2326 SETERRQ(1,"Must call MatStoreValues(A);first"); 2327 } 2328 2329 /* copy values over */ 2330 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2331 PetscFunctionReturn(0); 2332 } 2333 EXTERN_C_END 2334 2335 #undef __FUNCT__ 2336 #define __FUNCT__ "MatRetrieveValues" 2337 /*@ 2338 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2339 example, reuse of the linear part of a Jacobian, while recomputing the 2340 nonlinear portion. 2341 2342 Collect on Mat 2343 2344 Input Parameters: 2345 . mat - the matrix (currently on AIJ matrices support this option) 2346 2347 Level: advanced 2348 2349 .seealso: MatStoreValues() 2350 2351 @*/ 2352 int MatRetrieveValues(Mat mat) 2353 { 2354 int ierr,(*f)(Mat); 2355 2356 PetscFunctionBegin; 2357 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2358 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2359 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2360 2361 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2362 if (f) { 2363 ierr = (*f)(mat);CHKERRQ(ierr); 2364 } else { 2365 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2366 } 2367 PetscFunctionReturn(0); 2368 } 2369 2370 /* 2371 This allows SeqAIJ matrices to be passed to the matlab engine 2372 */ 2373 #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2374 #include "engine.h" /* Matlab include file */ 2375 #include "mex.h" /* Matlab include file */ 2376 EXTERN_C_BEGIN 2377 #undef __FUNCT__ 2378 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2379 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *mengine) 2380 { 2381 int ierr; 2382 Mat B = (Mat)obj; 2383 mxArray *mat; 2384 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2385 2386 PetscFunctionBegin; 2387 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2388 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2389 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2390 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2391 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2392 2393 /* Matlab indices start at 0 for sparse (what a surprise) */ 2394 2395 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2396 engPutVariable((Engine *)mengine,obj->name,mat); 2397 PetscFunctionReturn(0); 2398 } 2399 EXTERN_C_END 2400 2401 EXTERN_C_BEGIN 2402 #undef __FUNCT__ 2403 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2404 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *mengine) 2405 { 2406 int ierr,ii; 2407 Mat mat = (Mat)obj; 2408 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2409 mxArray *mmat; 2410 2411 PetscFunctionBegin; 2412 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2413 2414 mmat = engGetVariable((Engine *)mengine,obj->name); 2415 2416 aij->nz = (mxGetJc(mmat))[mat->m]; 2417 ierr = PetscMalloc(((size_t) aij->nz)*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2418 aij->j = (int*)(aij->a + aij->nz); 2419 aij->i = aij->j + aij->nz; 2420 aij->singlemalloc = PETSC_TRUE; 2421 aij->freedata = PETSC_TRUE; 2422 2423 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2424 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2425 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2426 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2427 2428 for (ii=0; ii<mat->m; ii++) { 2429 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2430 } 2431 2432 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2433 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2434 2435 PetscFunctionReturn(0); 2436 } 2437 EXTERN_C_END 2438 #endif 2439 2440 /* --------------------------------------------------------------------------------*/ 2441 #undef __FUNCT__ 2442 #define __FUNCT__ "MatCreateSeqAIJ" 2443 /*@C 2444 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2445 (the default parallel PETSc format). For good matrix assembly performance 2446 the user should preallocate the matrix storage by setting the parameter nz 2447 (or the array nnz). By setting these parameters accurately, performance 2448 during matrix assembly can be increased by more than a factor of 50. 2449 2450 Collective on MPI_Comm 2451 2452 Input Parameters: 2453 + comm - MPI communicator, set to PETSC_COMM_SELF 2454 . m - number of rows 2455 . n - number of columns 2456 . nz - number of nonzeros per row (same for all rows) 2457 - nnz - array containing the number of nonzeros in the various rows 2458 (possibly different for each row) or PETSC_NULL 2459 2460 Output Parameter: 2461 . A - the matrix 2462 2463 Notes: 2464 The AIJ format (also called the Yale sparse matrix format or 2465 compressed row storage), is fully compatible with standard Fortran 77 2466 storage. That is, the stored row and column indices can begin at 2467 either one (as in Fortran) or zero. See the users' manual for details. 2468 2469 Specify the preallocated storage with either nz or nnz (not both). 2470 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2471 allocation. For large problems you MUST preallocate memory or you 2472 will get TERRIBLE performance, see the users' manual chapter on matrices. 2473 2474 By default, this format uses inodes (identical nodes) when possible, to 2475 improve numerical efficiency of matrix-vector products and solves. We 2476 search for consecutive rows with the same nonzero structure, thereby 2477 reusing matrix information to achieve increased efficiency. 2478 2479 Options Database Keys: 2480 + -mat_aij_no_inode - Do not use inodes 2481 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2482 - -mat_aij_oneindex - Internally use indexing starting at 1 2483 rather than 0. Note that when calling MatSetValues(), 2484 the user still MUST index entries starting at 0! 2485 2486 Level: intermediate 2487 2488 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2489 2490 @*/ 2491 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,const int nnz[],Mat *A) 2492 { 2493 int ierr; 2494 2495 PetscFunctionBegin; 2496 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2497 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2498 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2499 PetscFunctionReturn(0); 2500 } 2501 2502 #define SKIP_ALLOCATION -4 2503 2504 #undef __FUNCT__ 2505 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2506 /*@C 2507 MatSeqAIJSetPreallocation - For good matrix assembly performance 2508 the user should preallocate the matrix storage by setting the parameter nz 2509 (or the array nnz). By setting these parameters accurately, performance 2510 during matrix assembly can be increased by more than a factor of 50. 2511 2512 Collective on MPI_Comm 2513 2514 Input Parameters: 2515 + comm - MPI communicator, set to PETSC_COMM_SELF 2516 . m - number of rows 2517 . n - number of columns 2518 . nz - number of nonzeros per row (same for all rows) 2519 - nnz - array containing the number of nonzeros in the various rows 2520 (possibly different for each row) or PETSC_NULL 2521 2522 Output Parameter: 2523 . A - the matrix 2524 2525 Notes: 2526 The AIJ format (also called the Yale sparse matrix format or 2527 compressed row storage), is fully compatible with standard Fortran 77 2528 storage. That is, the stored row and column indices can begin at 2529 either one (as in Fortran) or zero. See the users' manual for details. 2530 2531 Specify the preallocated storage with either nz or nnz (not both). 2532 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2533 allocation. For large problems you MUST preallocate memory or you 2534 will get TERRIBLE performance, see the users' manual chapter on matrices. 2535 2536 By default, this format uses inodes (identical nodes) when possible, to 2537 improve numerical efficiency of matrix-vector products and solves. We 2538 search for consecutive rows with the same nonzero structure, thereby 2539 reusing matrix information to achieve increased efficiency. 2540 2541 Options Database Keys: 2542 + -mat_aij_no_inode - Do not use inodes 2543 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2544 - -mat_aij_oneindex - Internally use indexing starting at 1 2545 rather than 0. Note that when calling MatSetValues(), 2546 the user still MUST index entries starting at 0! 2547 2548 Level: intermediate 2549 2550 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2551 2552 @*/ 2553 int MatSeqAIJSetPreallocation(Mat B,int nz,const int nnz[]) 2554 { 2555 int ierr,(*f)(Mat,int,const int[]); 2556 2557 PetscFunctionBegin; 2558 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2559 if (f) { 2560 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2561 } 2562 PetscFunctionReturn(0); 2563 } 2564 2565 EXTERN_C_BEGIN 2566 #undef __FUNCT__ 2567 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2568 int MatSeqAIJSetPreallocation_SeqAIJ(Mat B,int nz,int *nnz) 2569 { 2570 Mat_SeqAIJ *b; 2571 size_t len = 0; 2572 PetscTruth skipallocation = PETSC_FALSE; 2573 int i,ierr; 2574 2575 PetscFunctionBegin; 2576 2577 if (nz == SKIP_ALLOCATION) { 2578 skipallocation = PETSC_TRUE; 2579 nz = 0; 2580 } 2581 2582 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2583 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2584 if (nnz) { 2585 for (i=0; i<B->m; i++) { 2586 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2587 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); 2588 } 2589 } 2590 2591 B->preallocated = PETSC_TRUE; 2592 b = (Mat_SeqAIJ*)B->data; 2593 2594 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2595 if (!nnz) { 2596 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2597 else if (nz <= 0) nz = 1; 2598 for (i=0; i<B->m; i++) b->imax[i] = nz; 2599 nz = nz*B->m; 2600 } else { 2601 nz = 0; 2602 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2603 } 2604 2605 if (!skipallocation) { 2606 /* allocate the matrix space */ 2607 len = ((size_t) nz)*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2608 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2609 b->j = (int*)(b->a + nz); 2610 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2611 b->i = b->j + nz; 2612 b->i[0] = 0; 2613 for (i=1; i<B->m+1; i++) { 2614 b->i[i] = b->i[i-1] + b->imax[i-1]; 2615 } 2616 b->singlemalloc = PETSC_TRUE; 2617 b->freedata = PETSC_TRUE; 2618 } else { 2619 b->freedata = PETSC_FALSE; 2620 } 2621 2622 /* b->ilen will count nonzeros in each row so far. */ 2623 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2624 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2625 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2626 2627 b->nz = 0; 2628 b->maxnz = nz; 2629 B->info.nz_unneeded = (double)b->maxnz; 2630 PetscFunctionReturn(0); 2631 } 2632 EXTERN_C_END 2633 2634 /*MC 2635 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2636 based on compressed sparse row format. 2637 2638 Options Database Keys: 2639 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2640 2641 Level: beginner 2642 2643 .seealso: MatCreateSeqAIJ 2644 M*/ 2645 2646 EXTERN_C_BEGIN 2647 #undef __FUNCT__ 2648 #define __FUNCT__ "MatCreate_SeqAIJ" 2649 int MatCreate_SeqAIJ(Mat B) 2650 { 2651 Mat_SeqAIJ *b; 2652 int ierr,size; 2653 PetscTruth flg; 2654 2655 PetscFunctionBegin; 2656 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2657 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2658 2659 B->m = B->M = PetscMax(B->m,B->M); 2660 B->n = B->N = PetscMax(B->n,B->N); 2661 2662 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2663 B->data = (void*)b; 2664 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2665 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2666 B->factor = 0; 2667 B->lupivotthreshold = 1.0; 2668 B->mapping = 0; 2669 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2670 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2671 b->row = 0; 2672 b->col = 0; 2673 b->icol = 0; 2674 b->reallocs = 0; 2675 2676 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2677 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2678 2679 b->sorted = PETSC_FALSE; 2680 b->ignorezeroentries = PETSC_FALSE; 2681 b->roworiented = PETSC_TRUE; 2682 b->nonew = 0; 2683 b->diag = 0; 2684 b->solve_work = 0; 2685 B->spptr = 0; 2686 b->inode.use = PETSC_TRUE; 2687 b->inode.node_count = 0; 2688 b->inode.size = 0; 2689 b->inode.limit = 5; 2690 b->inode.max_limit = 5; 2691 b->saved_values = 0; 2692 b->idiag = 0; 2693 b->ssor = 0; 2694 b->keepzeroedrows = PETSC_FALSE; 2695 b->xtoy = 0; 2696 b->XtoY = 0; 2697 2698 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2699 2700 ierr = PetscOptionsHasName(B->prefix,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2701 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2702 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2703 "MatSeqAIJSetColumnIndices_SeqAIJ", 2704 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2705 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2706 "MatStoreValues_SeqAIJ", 2707 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2708 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2709 "MatRetrieveValues_SeqAIJ", 2710 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2711 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2712 "MatConvert_SeqAIJ_SeqSBAIJ", 2713 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2714 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2715 "MatConvert_SeqAIJ_SeqBAIJ", 2716 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2717 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2718 "MatIsTranspose_SeqAIJ", 2719 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2720 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2721 "MatSeqAIJSetPreallocation_SeqAIJ", 2722 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2723 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2724 "MatReorderForNonzeroDiagonal_SeqAIJ", 2725 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2726 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2727 "MatAdjustForInodes_SeqAIJ", 2728 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2729 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2730 "MatSeqAIJGetInodeSizes_SeqAIJ", 2731 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2732 #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2733 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2735 #endif 2736 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2737 PetscFunctionReturn(0); 2738 } 2739 EXTERN_C_END 2740 2741 #undef __FUNCT__ 2742 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2743 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2744 { 2745 Mat C; 2746 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2747 int i,m = A->m,ierr; 2748 size_t len; 2749 2750 PetscFunctionBegin; 2751 *B = 0; 2752 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2753 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2754 c = (Mat_SeqAIJ*)C->data; 2755 2756 C->factor = A->factor; 2757 c->row = 0; 2758 c->col = 0; 2759 c->icol = 0; 2760 c->keepzeroedrows = a->keepzeroedrows; 2761 C->assembled = PETSC_TRUE; 2762 2763 C->M = A->m; 2764 C->N = A->n; 2765 2766 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2767 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2768 for (i=0; i<m; i++) { 2769 c->imax[i] = a->imax[i]; 2770 c->ilen[i] = a->ilen[i]; 2771 } 2772 2773 /* allocate the matrix space */ 2774 c->singlemalloc = PETSC_TRUE; 2775 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2776 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2777 c->j = (int*)(c->a + a->i[m] ); 2778 c->i = c->j + a->i[m]; 2779 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2780 if (m > 0) { 2781 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(int));CHKERRQ(ierr); 2782 if (cpvalues == MAT_COPY_VALUES) { 2783 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2784 } else { 2785 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2786 } 2787 } 2788 2789 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2790 c->sorted = a->sorted; 2791 c->roworiented = a->roworiented; 2792 c->nonew = a->nonew; 2793 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2794 c->saved_values = 0; 2795 c->idiag = 0; 2796 c->ssor = 0; 2797 c->ignorezeroentries = a->ignorezeroentries; 2798 c->freedata = PETSC_TRUE; 2799 2800 if (a->diag) { 2801 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2802 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2803 for (i=0; i<m; i++) { 2804 c->diag[i] = a->diag[i]; 2805 } 2806 } else c->diag = 0; 2807 c->inode.use = a->inode.use; 2808 c->inode.limit = a->inode.limit; 2809 c->inode.max_limit = a->inode.max_limit; 2810 if (a->inode.size){ 2811 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2812 c->inode.node_count = a->inode.node_count; 2813 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2814 } else { 2815 c->inode.size = 0; 2816 c->inode.node_count = 0; 2817 } 2818 c->nz = a->nz; 2819 c->maxnz = a->maxnz; 2820 c->solve_work = 0; 2821 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2822 C->preallocated = PETSC_TRUE; 2823 2824 *B = C; 2825 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2826 PetscFunctionReturn(0); 2827 } 2828 2829 #undef __FUNCT__ 2830 #define __FUNCT__ "MatLoad_SeqAIJ" 2831 int MatLoad_SeqAIJ(PetscViewer viewer,const MatType type,Mat *A) 2832 { 2833 Mat_SeqAIJ *a; 2834 Mat B; 2835 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N; 2836 MPI_Comm comm; 2837 2838 PetscFunctionBegin; 2839 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2840 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2841 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2842 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2843 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2844 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2845 M = header[1]; N = header[2]; nz = header[3]; 2846 2847 if (nz < 0) { 2848 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2849 } 2850 2851 /* read in row lengths */ 2852 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2853 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2854 2855 /* create our matrix */ 2856 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2857 ierr = MatSetType(B,type);CHKERRQ(ierr); 2858 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2859 a = (Mat_SeqAIJ*)B->data; 2860 2861 /* read in column indices and adjust for Fortran indexing*/ 2862 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2863 2864 /* read in nonzero values */ 2865 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2866 2867 /* set matrix "i" values */ 2868 a->i[0] = 0; 2869 for (i=1; i<= M; i++) { 2870 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2871 a->ilen[i-1] = rowlengths[i-1]; 2872 } 2873 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2874 2875 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2876 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2877 *A = B; 2878 PetscFunctionReturn(0); 2879 } 2880 2881 #undef __FUNCT__ 2882 #define __FUNCT__ "MatEqual_SeqAIJ" 2883 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2884 { 2885 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2886 int ierr; 2887 2888 PetscFunctionBegin; 2889 /* If the matrix dimensions are not equal,or no of nonzeros */ 2890 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2891 *flg = PETSC_FALSE; 2892 PetscFunctionReturn(0); 2893 } 2894 2895 /* if the a->i are the same */ 2896 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2897 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2898 2899 /* if a->j are the same */ 2900 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2901 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2902 2903 /* if a->a are the same */ 2904 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2905 2906 PetscFunctionReturn(0); 2907 2908 } 2909 2910 #undef __FUNCT__ 2911 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2912 /*@C 2913 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2914 provided by the user. 2915 2916 Coolective on MPI_Comm 2917 2918 Input Parameters: 2919 + comm - must be an MPI communicator of size 1 2920 . m - number of rows 2921 . n - number of columns 2922 . i - row indices 2923 . j - column indices 2924 - a - matrix values 2925 2926 Output Parameter: 2927 . mat - the matrix 2928 2929 Level: intermediate 2930 2931 Notes: 2932 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2933 once the matrix is destroyed 2934 2935 You cannot set new nonzero locations into this matrix, that will generate an error. 2936 2937 The i and j indices are 0 based 2938 2939 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2940 2941 @*/ 2942 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2943 { 2944 int ierr,ii; 2945 Mat_SeqAIJ *aij; 2946 2947 PetscFunctionBegin; 2948 ierr = MatCreateSeqAIJ(comm,m,n,SKIP_ALLOCATION,0,mat);CHKERRQ(ierr); 2949 aij = (Mat_SeqAIJ*)(*mat)->data; 2950 2951 if (i[0] != 0) { 2952 SETERRQ(1,"i (row indices) must start with 0"); 2953 } 2954 aij->i = i; 2955 aij->j = j; 2956 aij->a = a; 2957 aij->singlemalloc = PETSC_FALSE; 2958 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2959 aij->freedata = PETSC_FALSE; 2960 2961 for (ii=0; ii<m; ii++) { 2962 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2963 #if defined(PETSC_USE_BOPT_g) 2964 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]); 2965 #endif 2966 } 2967 #if defined(PETSC_USE_BOPT_g) 2968 for (ii=0; ii<aij->i[m]; ii++) { 2969 if (j[ii] < 0) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2970 if (j[ii] > n - 1) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2971 } 2972 #endif 2973 2974 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2975 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2976 PetscFunctionReturn(0); 2977 } 2978 2979 #undef __FUNCT__ 2980 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2981 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2982 { 2983 int ierr; 2984 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2985 2986 PetscFunctionBegin; 2987 if (coloring->ctype == IS_COLORING_LOCAL) { 2988 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2989 a->coloring = coloring; 2990 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 2991 int i,*larray; 2992 ISColoring ocoloring; 2993 ISColoringValue *colors; 2994 2995 /* set coloring for diagonal portion */ 2996 ierr = PetscMalloc((A->n+1)*sizeof(int),&larray);CHKERRQ(ierr); 2997 for (i=0; i<A->n; i++) { 2998 larray[i] = i; 2999 } 3000 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3001 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3002 for (i=0; i<A->n; i++) { 3003 colors[i] = coloring->colors[larray[i]]; 3004 } 3005 ierr = PetscFree(larray);CHKERRQ(ierr); 3006 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3007 a->coloring = ocoloring; 3008 } 3009 PetscFunctionReturn(0); 3010 } 3011 3012 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 3013 EXTERN_C_BEGIN 3014 #include "adic/ad_utils.h" 3015 EXTERN_C_END 3016 3017 #undef __FUNCT__ 3018 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3019 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3020 { 3021 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3022 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3023 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3024 ISColoringValue *color; 3025 3026 PetscFunctionBegin; 3027 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3028 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3029 color = a->coloring->colors; 3030 /* loop over rows */ 3031 for (i=0; i<m; i++) { 3032 nz = ii[i+1] - ii[i]; 3033 /* loop over columns putting computed value into matrix */ 3034 for (j=0; j<nz; j++) { 3035 *v++ = values[color[*jj++]]; 3036 } 3037 values += nlen; /* jump to next row of derivatives */ 3038 } 3039 PetscFunctionReturn(0); 3040 } 3041 3042 #else 3043 3044 #undef __FUNCT__ 3045 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3046 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3047 { 3048 PetscFunctionBegin; 3049 SETERRQ(1,"PETSc installed without ADIC"); 3050 } 3051 3052 #endif 3053 3054 #undef __FUNCT__ 3055 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3056 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 3057 { 3058 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3059 int m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3060 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3061 ISColoringValue *color; 3062 3063 PetscFunctionBegin; 3064 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 3065 color = a->coloring->colors; 3066 /* loop over rows */ 3067 for (i=0; i<m; i++) { 3068 nz = ii[i+1] - ii[i]; 3069 /* loop over columns putting computed value into matrix */ 3070 for (j=0; j<nz; j++) { 3071 *v++ = values[color[*jj++]]; 3072 } 3073 values += nl; /* jump to next row of derivatives */ 3074 } 3075 PetscFunctionReturn(0); 3076 } 3077 3078