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 PetscFunctionReturn(0); 333 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 334 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 335 for (i=0; i<m; i++) { 336 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 337 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 338 #if defined(PETSC_USE_COMPLEX) 339 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 340 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 341 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 342 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 343 } else if (PetscRealPart(a->a[j]) != 0.0) { 344 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 345 } 346 #else 347 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 348 #endif 349 } 350 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 351 } 352 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 353 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 354 int nzd=0,fshift=1,*sptr; 355 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 356 ierr = PetscMalloc((m+1)*sizeof(int),&sptr);CHKERRQ(ierr); 357 for (i=0; i<m; i++) { 358 sptr[i] = nzd+1; 359 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 360 if (a->j[j] >= i) { 361 #if defined(PETSC_USE_COMPLEX) 362 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 363 #else 364 if (a->a[j] != 0.0) nzd++; 365 #endif 366 } 367 } 368 } 369 sptr[m] = nzd+1; 370 ierr = PetscViewerASCIIPrintf(viewer," %d %d\n\n",m,nzd);CHKERRQ(ierr); 371 for (i=0; i<m+1; i+=6) { 372 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);} 373 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);} 374 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);} 375 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 376 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 377 else {ierr = PetscViewerASCIIPrintf(viewer," %d\n",sptr[i]);CHKERRQ(ierr);} 378 } 379 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 380 ierr = PetscFree(sptr);CHKERRQ(ierr); 381 for (i=0; i<m; i++) { 382 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 383 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %d ",a->j[j]+fshift);CHKERRQ(ierr);} 384 } 385 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 386 } 387 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 388 for (i=0; i<m; i++) { 389 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 390 if (a->j[j] >= i) { 391 #if defined(PETSC_USE_COMPLEX) 392 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 393 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 394 } 395 #else 396 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 397 #endif 398 } 399 } 400 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 401 } 402 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 403 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 404 int cnt = 0,jcnt; 405 PetscScalar value; 406 407 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 408 for (i=0; i<m; i++) { 409 jcnt = 0; 410 for (j=0; j<A->n; j++) { 411 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 412 value = a->a[cnt++]; 413 jcnt++; 414 } else { 415 value = 0.0; 416 } 417 #if defined(PETSC_USE_COMPLEX) 418 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 419 #else 420 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 421 #endif 422 } 423 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 424 } 425 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 426 } else { 427 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 428 for (i=0; i<m; i++) { 429 ierr = PetscViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 430 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 431 #if defined(PETSC_USE_COMPLEX) 432 if (PetscImaginaryPart(a->a[j]) > 0.0) { 433 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 434 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 435 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 436 } else { 437 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 438 } 439 #else 440 ierr = PetscViewerASCIIPrintf(viewer," (%d, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 441 #endif 442 } 443 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 444 } 445 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 446 } 447 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 448 PetscFunctionReturn(0); 449 } 450 451 #undef __FUNCT__ 452 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 453 int MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 454 { 455 Mat A = (Mat) Aa; 456 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 457 int ierr,i,j,m = A->m,color; 458 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 459 PetscViewer viewer; 460 PetscViewerFormat format; 461 462 PetscFunctionBegin; 463 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 464 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 465 466 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 467 /* loop over matrix elements drawing boxes */ 468 469 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 470 /* Blue for negative, Cyan for zero and Red for positive */ 471 color = PETSC_DRAW_BLUE; 472 for (i=0; i<m; i++) { 473 y_l = m - i - 1.0; y_r = y_l + 1.0; 474 for (j=a->i[i]; j<a->i[i+1]; j++) { 475 x_l = a->j[j] ; x_r = x_l + 1.0; 476 #if defined(PETSC_USE_COMPLEX) 477 if (PetscRealPart(a->a[j]) >= 0.) continue; 478 #else 479 if (a->a[j] >= 0.) continue; 480 #endif 481 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 482 } 483 } 484 color = PETSC_DRAW_CYAN; 485 for (i=0; i<m; i++) { 486 y_l = m - i - 1.0; y_r = y_l + 1.0; 487 for (j=a->i[i]; j<a->i[i+1]; j++) { 488 x_l = a->j[j]; x_r = x_l + 1.0; 489 if (a->a[j] != 0.) continue; 490 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 491 } 492 } 493 color = PETSC_DRAW_RED; 494 for (i=0; i<m; i++) { 495 y_l = m - i - 1.0; y_r = y_l + 1.0; 496 for (j=a->i[i]; j<a->i[i+1]; j++) { 497 x_l = a->j[j]; x_r = x_l + 1.0; 498 #if defined(PETSC_USE_COMPLEX) 499 if (PetscRealPart(a->a[j]) <= 0.) continue; 500 #else 501 if (a->a[j] <= 0.) continue; 502 #endif 503 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 504 } 505 } 506 } else { 507 /* use contour shading to indicate magnitude of values */ 508 /* first determine max of all nonzero values */ 509 int nz = a->nz,count; 510 PetscDraw popup; 511 PetscReal scale; 512 513 for (i=0; i<nz; i++) { 514 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 515 } 516 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 517 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 518 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 519 count = 0; 520 for (i=0; i<m; i++) { 521 y_l = m - i - 1.0; y_r = y_l + 1.0; 522 for (j=a->i[i]; j<a->i[i+1]; j++) { 523 x_l = a->j[j]; x_r = x_l + 1.0; 524 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(a->a[count])); 525 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 526 count++; 527 } 528 } 529 } 530 PetscFunctionReturn(0); 531 } 532 533 #undef __FUNCT__ 534 #define __FUNCT__ "MatView_SeqAIJ_Draw" 535 int MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 536 { 537 int ierr; 538 PetscDraw draw; 539 PetscReal xr,yr,xl,yl,h,w; 540 PetscTruth isnull; 541 542 PetscFunctionBegin; 543 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 544 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 545 if (isnull) PetscFunctionReturn(0); 546 547 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 548 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 549 xr += w; yr += h; xl = -w; yl = -h; 550 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 551 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 552 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 556 #undef __FUNCT__ 557 #define __FUNCT__ "MatView_SeqAIJ" 558 int MatView_SeqAIJ(Mat A,PetscViewer viewer) 559 { 560 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 561 int ierr; 562 PetscTruth issocket,isascii,isbinary,isdraw; 563 564 PetscFunctionBegin; 565 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 566 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 567 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 568 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 569 if (issocket) { 570 ierr = PetscViewerSocketPutSparse_Private(viewer,A->m,A->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 571 } else if (isascii) { 572 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 573 } else if (isbinary) { 574 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 575 } else if (isdraw) { 576 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 577 } else { 578 SETERRQ1(1,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 579 } 580 PetscFunctionReturn(0); 581 } 582 583 #undef __FUNCT__ 584 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 585 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 586 { 587 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 588 int fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax,ierr; 589 int m = A->m,*ip,N,*ailen = a->ilen,rmax = 0; 590 PetscScalar *aa = a->a,*ap; 591 592 PetscFunctionBegin; 593 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 594 595 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 596 for (i=1; i<m; i++) { 597 /* move each row back by the amount of empty slots (fshift) before it*/ 598 fshift += imax[i-1] - ailen[i-1]; 599 rmax = PetscMax(rmax,ailen[i]); 600 if (fshift) { 601 ip = aj + ai[i] ; 602 ap = aa + ai[i] ; 603 N = ailen[i]; 604 for (j=0; j<N; j++) { 605 ip[j-fshift] = ip[j]; 606 ap[j-fshift] = ap[j]; 607 } 608 } 609 ai[i] = ai[i-1] + ailen[i-1]; 610 } 611 if (m) { 612 fshift += imax[m-1] - ailen[m-1]; 613 ai[m] = ai[m-1] + ailen[m-1]; 614 } 615 /* reset ilen and imax for each row */ 616 for (i=0; i<m; i++) { 617 ailen[i] = imax[i] = ai[i+1] - ai[i]; 618 } 619 a->nz = ai[m]; 620 621 /* diagonals may have moved, so kill the diagonal pointers */ 622 if (fshift && a->diag) { 623 ierr = PetscFree(a->diag);CHKERRQ(ierr); 624 PetscLogObjectMemory(A,-(m+1)*sizeof(int)); 625 a->diag = 0; 626 } 627 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded,%d used\n",m,A->n,fshift,a->nz); 628 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %d\n",a->reallocs); 629 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 630 a->reallocs = 0; 631 A->info.nz_unneeded = (double)fshift; 632 a->rmax = rmax; 633 634 /* check out for identical nodes. If found, use inode functions */ 635 ierr = Mat_AIJ_CheckInode(A,(PetscTruth)(!fshift));CHKERRQ(ierr); 636 637 PetscFunctionReturn(0); 638 } 639 640 #undef __FUNCT__ 641 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 642 int MatZeroEntries_SeqAIJ(Mat A) 643 { 644 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 645 int ierr; 646 647 PetscFunctionBegin; 648 ierr = PetscMemzero(a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 649 PetscFunctionReturn(0); 650 } 651 652 #undef __FUNCT__ 653 #define __FUNCT__ "MatDestroy_SeqAIJ" 654 int MatDestroy_SeqAIJ(Mat A) 655 { 656 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 657 int ierr; 658 659 PetscFunctionBegin; 660 #if defined(PETSC_USE_LOG) 661 PetscLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",A->m,A->n,a->nz); 662 #endif 663 if (a->freedata) { 664 ierr = PetscFree(a->a);CHKERRQ(ierr); 665 if (!a->singlemalloc) { 666 ierr = PetscFree(a->i);CHKERRQ(ierr); 667 ierr = PetscFree(a->j);CHKERRQ(ierr); 668 } 669 } 670 if (a->row) { 671 ierr = ISDestroy(a->row);CHKERRQ(ierr); 672 } 673 if (a->col) { 674 ierr = ISDestroy(a->col);CHKERRQ(ierr); 675 } 676 if (a->diag) {ierr = PetscFree(a->diag);CHKERRQ(ierr);} 677 if (a->ilen) {ierr = PetscFree(a->ilen);CHKERRQ(ierr);} 678 if (a->imax) {ierr = PetscFree(a->imax);CHKERRQ(ierr);} 679 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 680 if (a->solve_work) {ierr = PetscFree(a->solve_work);CHKERRQ(ierr);} 681 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 682 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 683 if (a->saved_values) {ierr = PetscFree(a->saved_values);CHKERRQ(ierr);} 684 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 685 if (a->xtoy) {ierr = PetscFree(a->xtoy);CHKERRQ(ierr);} 686 687 ierr = PetscFree(a);CHKERRQ(ierr); 688 PetscFunctionReturn(0); 689 } 690 691 #undef __FUNCT__ 692 #define __FUNCT__ "MatCompress_SeqAIJ" 693 int MatCompress_SeqAIJ(Mat A) 694 { 695 PetscFunctionBegin; 696 PetscFunctionReturn(0); 697 } 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "MatSetOption_SeqAIJ" 701 int MatSetOption_SeqAIJ(Mat A,MatOption op) 702 { 703 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 704 705 PetscFunctionBegin; 706 switch (op) { 707 case MAT_ROW_ORIENTED: 708 a->roworiented = PETSC_TRUE; 709 break; 710 case MAT_KEEP_ZEROED_ROWS: 711 a->keepzeroedrows = PETSC_TRUE; 712 break; 713 case MAT_COLUMN_ORIENTED: 714 a->roworiented = PETSC_FALSE; 715 break; 716 case MAT_COLUMNS_SORTED: 717 a->sorted = PETSC_TRUE; 718 break; 719 case MAT_COLUMNS_UNSORTED: 720 a->sorted = PETSC_FALSE; 721 break; 722 case MAT_NO_NEW_NONZERO_LOCATIONS: 723 a->nonew = 1; 724 break; 725 case MAT_NEW_NONZERO_LOCATION_ERR: 726 a->nonew = -1; 727 break; 728 case MAT_NEW_NONZERO_ALLOCATION_ERR: 729 a->nonew = -2; 730 break; 731 case MAT_YES_NEW_NONZERO_LOCATIONS: 732 a->nonew = 0; 733 break; 734 case MAT_IGNORE_ZERO_ENTRIES: 735 a->ignorezeroentries = PETSC_TRUE; 736 break; 737 case MAT_USE_INODES: 738 a->inode.use = PETSC_TRUE; 739 break; 740 case MAT_DO_NOT_USE_INODES: 741 a->inode.use = PETSC_FALSE; 742 break; 743 case MAT_ROWS_SORTED: 744 case MAT_ROWS_UNSORTED: 745 case MAT_YES_NEW_DIAGONALS: 746 case MAT_IGNORE_OFF_PROC_ENTRIES: 747 case MAT_USE_HASH_TABLE: 748 PetscLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 749 break; 750 case MAT_NO_NEW_DIAGONALS: 751 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 752 case MAT_INODE_LIMIT_1: 753 a->inode.limit = 1; 754 break; 755 case MAT_INODE_LIMIT_2: 756 a->inode.limit = 2; 757 break; 758 case MAT_INODE_LIMIT_3: 759 a->inode.limit = 3; 760 break; 761 case MAT_INODE_LIMIT_4: 762 a->inode.limit = 4; 763 break; 764 case MAT_INODE_LIMIT_5: 765 a->inode.limit = 5; 766 break; 767 case MAT_SYMMETRIC: 768 case MAT_STRUCTURALLY_SYMMETRIC: 769 case MAT_NOT_SYMMETRIC: 770 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 771 case MAT_HERMITIAN: 772 case MAT_NOT_HERMITIAN: 773 case MAT_SYMMETRY_ETERNAL: 774 case MAT_NOT_SYMMETRY_ETERNAL: 775 break; 776 default: 777 SETERRQ(PETSC_ERR_SUP,"unknown option"); 778 } 779 PetscFunctionReturn(0); 780 } 781 782 #undef __FUNCT__ 783 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 784 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 785 { 786 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 787 int i,j,n,ierr; 788 PetscScalar *x,zero = 0.0; 789 790 PetscFunctionBegin; 791 ierr = VecSet(&zero,v);CHKERRQ(ierr); 792 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 793 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 794 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 795 for (i=0; i<A->m; i++) { 796 for (j=a->i[i]; j<a->i[i+1]; j++) { 797 if (a->j[j] == i) { 798 x[i] = a->a[j]; 799 break; 800 } 801 } 802 } 803 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 804 PetscFunctionReturn(0); 805 } 806 807 808 #undef __FUNCT__ 809 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 810 int MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 811 { 812 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 813 PetscScalar *x,*y; 814 int ierr,m = A->m; 815 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 816 PetscScalar *v,alpha; 817 int n,i,*idx; 818 #endif 819 820 PetscFunctionBegin; 821 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 822 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 823 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 824 825 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 826 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 827 #else 828 for (i=0; i<m; i++) { 829 idx = a->j + a->i[i] ; 830 v = a->a + a->i[i] ; 831 n = a->i[i+1] - a->i[i]; 832 alpha = x[i]; 833 while (n-->0) {y[*idx++] += alpha * *v++;} 834 } 835 #endif 836 PetscLogFlops(2*a->nz); 837 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 838 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 839 PetscFunctionReturn(0); 840 } 841 842 #undef __FUNCT__ 843 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 844 int MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 845 { 846 PetscScalar zero = 0.0; 847 int ierr; 848 849 PetscFunctionBegin; 850 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 851 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 852 PetscFunctionReturn(0); 853 } 854 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "MatMult_SeqAIJ" 858 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 859 { 860 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 861 PetscScalar *x,*y,*v; 862 int ierr,m = A->m,*idx,*ii; 863 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 864 int n,i,jrow,j; 865 PetscScalar sum; 866 #endif 867 868 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 869 #pragma disjoint(*x,*y,*v) 870 #endif 871 872 PetscFunctionBegin; 873 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 874 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 875 idx = a->j; 876 v = a->a; 877 ii = a->i; 878 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 879 fortranmultaij_(&m,x,ii,idx,v,y); 880 #else 881 for (i=0; i<m; i++) { 882 jrow = ii[i]; 883 n = ii[i+1] - jrow; 884 sum = 0.0; 885 for (j=0; j<n; j++) { 886 sum += v[jrow]*x[idx[jrow]]; jrow++; 887 } 888 y[i] = sum; 889 } 890 #endif 891 PetscLogFlops(2*a->nz - m); 892 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 893 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 894 PetscFunctionReturn(0); 895 } 896 897 #undef __FUNCT__ 898 #define __FUNCT__ "MatMultAdd_SeqAIJ" 899 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 900 { 901 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 902 PetscScalar *x,*y,*z,*v; 903 int ierr,m = A->m,*idx,*ii; 904 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 905 int n,i,jrow,j; 906 PetscScalar sum; 907 #endif 908 909 PetscFunctionBegin; 910 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 911 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 912 if (zz != yy) { 913 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 914 } else { 915 z = y; 916 } 917 918 idx = a->j; 919 v = a->a; 920 ii = a->i; 921 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 922 fortranmultaddaij_(&m,x,ii,idx,v,y,z); 923 #else 924 for (i=0; i<m; i++) { 925 jrow = ii[i]; 926 n = ii[i+1] - jrow; 927 sum = y[i]; 928 for (j=0; j<n; j++) { 929 sum += v[jrow]*x[idx[jrow]]; jrow++; 930 } 931 z[i] = sum; 932 } 933 #endif 934 PetscLogFlops(2*a->nz); 935 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 936 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 937 if (zz != yy) { 938 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 939 } 940 PetscFunctionReturn(0); 941 } 942 943 /* 944 Adds diagonal pointers to sparse matrix structure. 945 */ 946 #undef __FUNCT__ 947 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 948 int MatMarkDiagonal_SeqAIJ(Mat A) 949 { 950 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 951 int i,j,*diag,m = A->m,ierr; 952 953 PetscFunctionBegin; 954 if (a->diag) PetscFunctionReturn(0); 955 956 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 957 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 958 for (i=0; i<A->m; i++) { 959 diag[i] = a->i[i+1]; 960 for (j=a->i[i]; j<a->i[i+1]; j++) { 961 if (a->j[j] == i) { 962 diag[i] = j; 963 break; 964 } 965 } 966 } 967 a->diag = diag; 968 PetscFunctionReturn(0); 969 } 970 971 /* 972 Checks for missing diagonals 973 */ 974 #undef __FUNCT__ 975 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 976 int MatMissingDiagonal_SeqAIJ(Mat A) 977 { 978 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 979 int *diag,*jj = a->j,i,ierr; 980 981 PetscFunctionBegin; 982 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 983 diag = a->diag; 984 for (i=0; i<A->m; i++) { 985 if (jj[diag[i]] != i) { 986 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 987 } 988 } 989 PetscFunctionReturn(0); 990 } 991 992 #undef __FUNCT__ 993 #define __FUNCT__ "MatRelax_SeqAIJ" 994 int MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,int lits,Vec xx) 995 { 996 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 997 PetscScalar *x,d,*xs,sum,*t,scale,*idiag=0,*mdiag; 998 const PetscScalar *v = a->a, *b, *bs,*xb, *ts; 999 int ierr,n = A->n,m = A->m,i; 1000 const int *idx,*diag; 1001 1002 PetscFunctionBegin; 1003 its = its*lits; 1004 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %d and local its %d both positive",its,lits); 1005 1006 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1007 diag = a->diag; 1008 if (!a->idiag) { 1009 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1010 a->ssor = a->idiag + m; 1011 mdiag = a->ssor + m; 1012 1013 v = a->a; 1014 1015 /* this is wrong when fshift omega changes each iteration */ 1016 if (omega == 1.0 && fshift == 0.0) { 1017 for (i=0; i<m; i++) { 1018 mdiag[i] = v[diag[i]]; 1019 a->idiag[i] = 1.0/v[diag[i]]; 1020 } 1021 PetscLogFlops(m); 1022 } else { 1023 for (i=0; i<m; i++) { 1024 mdiag[i] = v[diag[i]]; 1025 a->idiag[i] = omega/(fshift + v[diag[i]]); 1026 } 1027 PetscLogFlops(2*m); 1028 } 1029 } 1030 t = a->ssor; 1031 idiag = a->idiag; 1032 mdiag = a->idiag + 2*m; 1033 1034 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1035 if (xx != bb) { 1036 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1037 } else { 1038 b = x; 1039 } 1040 1041 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1042 xs = x; 1043 if (flag == SOR_APPLY_UPPER) { 1044 /* apply (U + D/omega) to the vector */ 1045 bs = b; 1046 for (i=0; i<m; i++) { 1047 d = fshift + a->a[diag[i]]; 1048 n = a->i[i+1] - diag[i] - 1; 1049 idx = a->j + diag[i] + 1; 1050 v = a->a + diag[i] + 1; 1051 sum = b[i]*d/omega; 1052 SPARSEDENSEDOT(sum,bs,v,idx,n); 1053 x[i] = sum; 1054 } 1055 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1056 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1057 PetscLogFlops(a->nz); 1058 PetscFunctionReturn(0); 1059 } 1060 1061 1062 /* Let A = L + U + D; where L is lower trianglar, 1063 U is upper triangular, E is diagonal; This routine applies 1064 1065 (L + E)^{-1} A (U + E)^{-1} 1066 1067 to a vector efficiently using Eisenstat's trick. This is for 1068 the case of SSOR preconditioner, so E is D/omega where omega 1069 is the relaxation factor. 1070 */ 1071 1072 if (flag == SOR_APPLY_LOWER) { 1073 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1074 } else if (flag & SOR_EISENSTAT) { 1075 /* Let A = L + U + D; where L is lower trianglar, 1076 U is upper triangular, E is diagonal; This routine applies 1077 1078 (L + E)^{-1} A (U + E)^{-1} 1079 1080 to a vector efficiently using Eisenstat's trick. This is for 1081 the case of SSOR preconditioner, so E is D/omega where omega 1082 is the relaxation factor. 1083 */ 1084 scale = (2.0/omega) - 1.0; 1085 1086 /* x = (E + U)^{-1} b */ 1087 for (i=m-1; i>=0; i--) { 1088 n = a->i[i+1] - diag[i] - 1; 1089 idx = a->j + diag[i] + 1; 1090 v = a->a + diag[i] + 1; 1091 sum = b[i]; 1092 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1093 x[i] = sum*idiag[i]; 1094 } 1095 1096 /* t = b - (2*E - D)x */ 1097 v = a->a; 1098 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1099 1100 /* t = (E + L)^{-1}t */ 1101 ts = t; 1102 diag = a->diag; 1103 for (i=0; i<m; i++) { 1104 n = diag[i] - a->i[i]; 1105 idx = a->j + a->i[i]; 1106 v = a->a + a->i[i]; 1107 sum = t[i]; 1108 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1109 t[i] = sum*idiag[i]; 1110 /* x = x + t */ 1111 x[i] += t[i]; 1112 } 1113 1114 PetscLogFlops(6*m-1 + 2*a->nz); 1115 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1116 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1117 PetscFunctionReturn(0); 1118 } 1119 if (flag & SOR_ZERO_INITIAL_GUESS) { 1120 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1121 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1122 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(int *)diag,idiag,a->a,(void*)b); 1123 #else 1124 for (i=0; i<m; i++) { 1125 n = diag[i] - a->i[i]; 1126 idx = a->j + a->i[i]; 1127 v = a->a + a->i[i]; 1128 sum = b[i]; 1129 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1130 x[i] = sum*idiag[i]; 1131 } 1132 #endif 1133 xb = x; 1134 PetscLogFlops(a->nz); 1135 } else xb = b; 1136 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1137 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1138 for (i=0; i<m; i++) { 1139 x[i] *= mdiag[i]; 1140 } 1141 PetscLogFlops(m); 1142 } 1143 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1144 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1145 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(int*)diag,idiag,a->a,(void*)xb); 1146 #else 1147 for (i=m-1; i>=0; i--) { 1148 n = a->i[i+1] - diag[i] - 1; 1149 idx = a->j + diag[i] + 1; 1150 v = a->a + diag[i] + 1; 1151 sum = xb[i]; 1152 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1153 x[i] = sum*idiag[i]; 1154 } 1155 #endif 1156 PetscLogFlops(a->nz); 1157 } 1158 its--; 1159 } 1160 while (its--) { 1161 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1162 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1163 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(int*)diag,a->a,(void*)b); 1164 #else 1165 for (i=0; i<m; i++) { 1166 d = fshift + a->a[diag[i]]; 1167 n = a->i[i+1] - a->i[i]; 1168 idx = a->j + a->i[i]; 1169 v = a->a + a->i[i]; 1170 sum = b[i]; 1171 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1172 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1173 } 1174 #endif 1175 PetscLogFlops(a->nz); 1176 } 1177 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1178 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1179 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(int*)diag,a->a,(void*)b); 1180 #else 1181 for (i=m-1; i>=0; i--) { 1182 d = fshift + a->a[diag[i]]; 1183 n = a->i[i+1] - a->i[i]; 1184 idx = a->j + a->i[i]; 1185 v = a->a + a->i[i]; 1186 sum = b[i]; 1187 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1188 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1189 } 1190 #endif 1191 PetscLogFlops(a->nz); 1192 } 1193 } 1194 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1195 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1196 PetscFunctionReturn(0); 1197 } 1198 1199 #undef __FUNCT__ 1200 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1201 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1202 { 1203 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1204 1205 PetscFunctionBegin; 1206 info->rows_global = (double)A->m; 1207 info->columns_global = (double)A->n; 1208 info->rows_local = (double)A->m; 1209 info->columns_local = (double)A->n; 1210 info->block_size = 1.0; 1211 info->nz_allocated = (double)a->maxnz; 1212 info->nz_used = (double)a->nz; 1213 info->nz_unneeded = (double)(a->maxnz - a->nz); 1214 info->assemblies = (double)A->num_ass; 1215 info->mallocs = (double)a->reallocs; 1216 info->memory = A->mem; 1217 if (A->factor) { 1218 info->fill_ratio_given = A->info.fill_ratio_given; 1219 info->fill_ratio_needed = A->info.fill_ratio_needed; 1220 info->factor_mallocs = A->info.factor_mallocs; 1221 } else { 1222 info->fill_ratio_given = 0; 1223 info->fill_ratio_needed = 0; 1224 info->factor_mallocs = 0; 1225 } 1226 PetscFunctionReturn(0); 1227 } 1228 1229 #undef __FUNCT__ 1230 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1231 int MatZeroRows_SeqAIJ(Mat A,IS is,const PetscScalar *diag) 1232 { 1233 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1234 int i,ierr,N,*rows,m = A->m - 1; 1235 1236 PetscFunctionBegin; 1237 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1238 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1239 if (a->keepzeroedrows) { 1240 for (i=0; i<N; i++) { 1241 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %d out of range", rows[i]); 1242 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1243 } 1244 if (diag) { 1245 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1246 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1247 for (i=0; i<N; i++) { 1248 a->a[a->diag[rows[i]]] = *diag; 1249 } 1250 } 1251 } else { 1252 if (diag) { 1253 for (i=0; i<N; i++) { 1254 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %d out of range", rows[i]); 1255 if (a->ilen[rows[i]] > 0) { 1256 a->ilen[rows[i]] = 1; 1257 a->a[a->i[rows[i]]] = *diag; 1258 a->j[a->i[rows[i]]] = rows[i]; 1259 } else { /* in case row was completely empty */ 1260 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1261 } 1262 } 1263 } else { 1264 for (i=0; i<N; i++) { 1265 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %d out of range", rows[i]); 1266 a->ilen[rows[i]] = 0; 1267 } 1268 } 1269 } 1270 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1271 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1272 PetscFunctionReturn(0); 1273 } 1274 1275 #undef __FUNCT__ 1276 #define __FUNCT__ "MatGetRow_SeqAIJ" 1277 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1278 { 1279 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1280 int *itmp; 1281 1282 PetscFunctionBegin; 1283 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %d out of range",row); 1284 1285 *nz = a->i[row+1] - a->i[row]; 1286 if (v) *v = a->a + a->i[row]; 1287 if (idx) { 1288 itmp = a->j + a->i[row]; 1289 if (*nz) { 1290 *idx = itmp; 1291 } 1292 else *idx = 0; 1293 } 1294 PetscFunctionReturn(0); 1295 } 1296 1297 /* remove this function? */ 1298 #undef __FUNCT__ 1299 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1300 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1301 { 1302 /* Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1303 int ierr; */ 1304 1305 PetscFunctionBegin; 1306 /* if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} */ 1307 PetscFunctionReturn(0); 1308 } 1309 1310 #undef __FUNCT__ 1311 #define __FUNCT__ "MatNorm_SeqAIJ" 1312 int MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1313 { 1314 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1315 PetscScalar *v = a->a; 1316 PetscReal sum = 0.0; 1317 int i,j,ierr; 1318 1319 PetscFunctionBegin; 1320 if (type == NORM_FROBENIUS) { 1321 for (i=0; i<a->nz; i++) { 1322 #if defined(PETSC_USE_COMPLEX) 1323 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1324 #else 1325 sum += (*v)*(*v); v++; 1326 #endif 1327 } 1328 *nrm = sqrt(sum); 1329 } else if (type == NORM_1) { 1330 PetscReal *tmp; 1331 int *jj = a->j; 1332 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1333 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1334 *nrm = 0.0; 1335 for (j=0; j<a->nz; j++) { 1336 tmp[*jj++] += PetscAbsScalar(*v); v++; 1337 } 1338 for (j=0; j<A->n; j++) { 1339 if (tmp[j] > *nrm) *nrm = tmp[j]; 1340 } 1341 ierr = PetscFree(tmp);CHKERRQ(ierr); 1342 } else if (type == NORM_INFINITY) { 1343 *nrm = 0.0; 1344 for (j=0; j<A->m; j++) { 1345 v = a->a + a->i[j]; 1346 sum = 0.0; 1347 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1348 sum += PetscAbsScalar(*v); v++; 1349 } 1350 if (sum > *nrm) *nrm = sum; 1351 } 1352 } else { 1353 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1354 } 1355 PetscFunctionReturn(0); 1356 } 1357 1358 #undef __FUNCT__ 1359 #define __FUNCT__ "MatTranspose_SeqAIJ" 1360 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1361 { 1362 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1363 Mat C; 1364 int i,ierr,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1365 PetscScalar *array = a->a; 1366 1367 PetscFunctionBegin; 1368 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1369 ierr = PetscMalloc((1+A->n)*sizeof(int),&col);CHKERRQ(ierr); 1370 ierr = PetscMemzero(col,(1+A->n)*sizeof(int));CHKERRQ(ierr); 1371 1372 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1373 ierr = MatCreate(A->comm,A->n,m,A->n,m,&C);CHKERRQ(ierr); 1374 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1375 ierr = MatSeqAIJSetPreallocation(C,0,col);CHKERRQ(ierr); 1376 ierr = PetscFree(col);CHKERRQ(ierr); 1377 for (i=0; i<m; i++) { 1378 len = ai[i+1]-ai[i]; 1379 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1380 array += len; 1381 aj += len; 1382 } 1383 1384 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1385 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1386 1387 if (B) { 1388 *B = C; 1389 } else { 1390 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1391 } 1392 PetscFunctionReturn(0); 1393 } 1394 1395 EXTERN_C_BEGIN 1396 #undef __FUNCT__ 1397 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1398 int MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscTruth *f) 1399 { 1400 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1401 int *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1402 int ma,na,mb,nb, i,ierr; 1403 1404 PetscFunctionBegin; 1405 bij = (Mat_SeqAIJ *) B->data; 1406 1407 ierr = MatGetSize(A,&ma,&na); CHKERRQ(ierr); 1408 ierr = MatGetSize(B,&mb,&nb); CHKERRQ(ierr); 1409 if (ma!=nb || na!=mb) 1410 SETERRQ(1,"Incompatible A/B sizes for symmetry test"); 1411 aii = aij->i; bii = bij->i; 1412 adx = aij->j; bdx = bij->j; 1413 va = aij->a; vb = bij->a; 1414 ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr); 1415 ierr = PetscMalloc(mb*sizeof(int),&bptr); CHKERRQ(ierr); 1416 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1417 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1418 1419 *f = PETSC_TRUE; 1420 for (i=0; i<ma; i++) { 1421 /*printf("row %d spans %d--%d; we start @ %d\n", 1422 i,idx[ii[i]],idx[ii[i+1]-1],idx[aptr[i]]);*/ 1423 while (aptr[i]<aii[i+1]) { 1424 int idc,idr; PetscScalar vc,vr; 1425 /* column/row index/value */ 1426 idc = adx[aptr[i]]; idr = bdx[bptr[idc]]; 1427 vc = va[aptr[i]]; vr = vb[bptr[idc]]; 1428 /*printf("comparing %d: (%d,%d)=%e to (%d,%d)=%e\n", 1429 aptr[i],i,idc,vc,idc,idr,vr);*/ 1430 if (i!=idr || vc!=vr) { 1431 *f = PETSC_FALSE; goto done; 1432 } else { 1433 aptr[i]++; if (B || i!=idc) bptr[idc]++; 1434 } 1435 } 1436 } 1437 done: 1438 ierr = PetscFree(aptr); CHKERRQ(ierr); 1439 if (B) { 1440 ierr = PetscFree(bptr); CHKERRQ(ierr); 1441 } 1442 1443 PetscFunctionReturn(0); 1444 } 1445 EXTERN_C_END 1446 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1449 int MatIsSymmetric_SeqAIJ(Mat A,PetscTruth *f) 1450 { 1451 int ierr; 1452 PetscFunctionBegin; 1453 ierr = MatIsTranspose_SeqAIJ(A,A,f); CHKERRQ(ierr); 1454 PetscFunctionReturn(0); 1455 } 1456 1457 #undef __FUNCT__ 1458 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1459 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1460 { 1461 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1462 PetscScalar *l,*r,x,*v; 1463 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj; 1464 1465 PetscFunctionBegin; 1466 if (ll) { 1467 /* The local size is used so that VecMPI can be passed to this routine 1468 by MatDiagonalScale_MPIAIJ */ 1469 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1470 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1471 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1472 v = a->a; 1473 for (i=0; i<m; i++) { 1474 x = l[i]; 1475 M = a->i[i+1] - a->i[i]; 1476 for (j=0; j<M; j++) { (*v++) *= x;} 1477 } 1478 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1479 PetscLogFlops(nz); 1480 } 1481 if (rr) { 1482 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1483 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1484 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1485 v = a->a; jj = a->j; 1486 for (i=0; i<nz; i++) { 1487 (*v++) *= r[*jj++]; 1488 } 1489 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1490 PetscLogFlops(nz); 1491 } 1492 PetscFunctionReturn(0); 1493 } 1494 1495 #undef __FUNCT__ 1496 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1497 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1498 { 1499 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1500 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1501 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1502 int *irow,*icol,nrows,ncols; 1503 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1504 PetscScalar *a_new,*mat_a; 1505 Mat C; 1506 PetscTruth stride; 1507 1508 PetscFunctionBegin; 1509 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1510 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1511 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1512 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1513 1514 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1515 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1516 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1517 1518 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1519 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1520 if (stride && step == 1) { 1521 /* special case of contiguous rows */ 1522 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1523 starts = lens + nrows; 1524 /* loop over new rows determining lens and starting points */ 1525 for (i=0; i<nrows; i++) { 1526 kstart = ai[irow[i]]; 1527 kend = kstart + ailen[irow[i]]; 1528 for (k=kstart; k<kend; k++) { 1529 if (aj[k] >= first) { 1530 starts[i] = k; 1531 break; 1532 } 1533 } 1534 sum = 0; 1535 while (k < kend) { 1536 if (aj[k++] >= first+ncols) break; 1537 sum++; 1538 } 1539 lens[i] = sum; 1540 } 1541 /* create submatrix */ 1542 if (scall == MAT_REUSE_MATRIX) { 1543 int n_cols,n_rows; 1544 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1545 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1546 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1547 C = *B; 1548 } else { 1549 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1550 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1551 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1552 } 1553 c = (Mat_SeqAIJ*)C->data; 1554 1555 /* loop over rows inserting into submatrix */ 1556 a_new = c->a; 1557 j_new = c->j; 1558 i_new = c->i; 1559 1560 for (i=0; i<nrows; i++) { 1561 ii = starts[i]; 1562 lensi = lens[i]; 1563 for (k=0; k<lensi; k++) { 1564 *j_new++ = aj[ii+k] - first; 1565 } 1566 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1567 a_new += lensi; 1568 i_new[i+1] = i_new[i] + lensi; 1569 c->ilen[i] = lensi; 1570 } 1571 ierr = PetscFree(lens);CHKERRQ(ierr); 1572 } else { 1573 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1574 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1575 1576 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1577 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1578 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1579 /* determine lens of each row */ 1580 for (i=0; i<nrows; i++) { 1581 kstart = ai[irow[i]]; 1582 kend = kstart + a->ilen[irow[i]]; 1583 lens[i] = 0; 1584 for (k=kstart; k<kend; k++) { 1585 if (smap[aj[k]]) { 1586 lens[i]++; 1587 } 1588 } 1589 } 1590 /* Create and fill new matrix */ 1591 if (scall == MAT_REUSE_MATRIX) { 1592 PetscTruth equal; 1593 1594 c = (Mat_SeqAIJ *)((*B)->data); 1595 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1596 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1597 if (!equal) { 1598 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1599 } 1600 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1601 C = *B; 1602 } else { 1603 ierr = MatCreate(A->comm,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE,&C);CHKERRQ(ierr); 1604 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1605 ierr = MatSeqAIJSetPreallocation(C,0,lens);CHKERRQ(ierr); 1606 } 1607 c = (Mat_SeqAIJ *)(C->data); 1608 for (i=0; i<nrows; i++) { 1609 row = irow[i]; 1610 kstart = ai[row]; 1611 kend = kstart + a->ilen[row]; 1612 mat_i = c->i[i]; 1613 mat_j = c->j + mat_i; 1614 mat_a = c->a + mat_i; 1615 mat_ilen = c->ilen + i; 1616 for (k=kstart; k<kend; k++) { 1617 if ((tcol=smap[a->j[k]])) { 1618 *mat_j++ = tcol - 1; 1619 *mat_a++ = a->a[k]; 1620 (*mat_ilen)++; 1621 1622 } 1623 } 1624 } 1625 /* Free work space */ 1626 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1627 ierr = PetscFree(smap);CHKERRQ(ierr); 1628 ierr = PetscFree(lens);CHKERRQ(ierr); 1629 } 1630 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1631 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1632 1633 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1634 *B = C; 1635 PetscFunctionReturn(0); 1636 } 1637 1638 /* 1639 */ 1640 #undef __FUNCT__ 1641 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1642 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1643 { 1644 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1645 int ierr; 1646 Mat outA; 1647 PetscTruth row_identity,col_identity; 1648 1649 PetscFunctionBegin; 1650 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1651 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1652 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1653 if (!row_identity || !col_identity) { 1654 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1655 } 1656 1657 outA = inA; 1658 inA->factor = FACTOR_LU; 1659 a->row = row; 1660 a->col = col; 1661 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1662 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1663 1664 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1665 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1666 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1667 PetscLogObjectParent(inA,a->icol); 1668 1669 if (!a->solve_work) { /* this matrix may have been factored before */ 1670 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1671 } 1672 1673 if (!a->diag) { 1674 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1675 } 1676 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1677 PetscFunctionReturn(0); 1678 } 1679 1680 #include "petscblaslapack.h" 1681 #undef __FUNCT__ 1682 #define __FUNCT__ "MatScale_SeqAIJ" 1683 int MatScale_SeqAIJ(const PetscScalar *alpha,Mat inA) 1684 { 1685 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1686 int one = 1; 1687 1688 PetscFunctionBegin; 1689 BLscal_(&a->nz,(PetscScalar*)alpha,a->a,&one); 1690 PetscLogFlops(a->nz); 1691 PetscFunctionReturn(0); 1692 } 1693 1694 #undef __FUNCT__ 1695 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1696 int MatGetSubMatrices_SeqAIJ(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1697 { 1698 int ierr,i; 1699 1700 PetscFunctionBegin; 1701 if (scall == MAT_INITIAL_MATRIX) { 1702 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1703 } 1704 1705 for (i=0; i<n; i++) { 1706 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1707 } 1708 PetscFunctionReturn(0); 1709 } 1710 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1713 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1714 { 1715 PetscFunctionBegin; 1716 *bs = 1; 1717 PetscFunctionReturn(0); 1718 } 1719 1720 #undef __FUNCT__ 1721 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1722 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS is[],int ov) 1723 { 1724 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1725 int row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1726 int start,end,*ai,*aj; 1727 PetscBT table; 1728 1729 PetscFunctionBegin; 1730 m = A->m; 1731 ai = a->i; 1732 aj = a->j; 1733 1734 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1735 1736 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1737 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1738 1739 for (i=0; i<is_max; i++) { 1740 /* Initialize the two local arrays */ 1741 isz = 0; 1742 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1743 1744 /* Extract the indices, assume there can be duplicate entries */ 1745 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1746 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1747 1748 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1749 for (j=0; j<n ; ++j){ 1750 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1751 } 1752 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1753 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1754 1755 k = 0; 1756 for (j=0; j<ov; j++){ /* for each overlap */ 1757 n = isz; 1758 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1759 row = nidx[k]; 1760 start = ai[row]; 1761 end = ai[row+1]; 1762 for (l = start; l<end ; l++){ 1763 val = aj[l] ; 1764 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1765 } 1766 } 1767 } 1768 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1769 } 1770 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1771 ierr = PetscFree(nidx);CHKERRQ(ierr); 1772 PetscFunctionReturn(0); 1773 } 1774 1775 /* -------------------------------------------------------------- */ 1776 #undef __FUNCT__ 1777 #define __FUNCT__ "MatPermute_SeqAIJ" 1778 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1779 { 1780 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1781 PetscScalar *vwork; 1782 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1783 int *row,*col,*cnew,j,*lens; 1784 IS icolp,irowp; 1785 1786 PetscFunctionBegin; 1787 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1788 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1789 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1790 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1791 1792 /* determine lengths of permuted rows */ 1793 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1794 for (i=0; i<m; i++) { 1795 lens[row[i]] = a->i[i+1] - a->i[i]; 1796 } 1797 ierr = MatCreate(A->comm,m,n,m,n,B);CHKERRQ(ierr); 1798 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 1799 ierr = MatSeqAIJSetPreallocation(*B,0,lens);CHKERRQ(ierr); 1800 ierr = PetscFree(lens);CHKERRQ(ierr); 1801 1802 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1803 for (i=0; i<m; i++) { 1804 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1805 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1806 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1807 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1808 } 1809 ierr = PetscFree(cnew);CHKERRQ(ierr); 1810 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1811 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1812 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1813 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1814 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1815 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1816 PetscFunctionReturn(0); 1817 } 1818 1819 #undef __FUNCT__ 1820 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1821 int MatPrintHelp_SeqAIJ(Mat A) 1822 { 1823 static PetscTruth called = PETSC_FALSE; 1824 MPI_Comm comm = A->comm; 1825 int ierr; 1826 1827 PetscFunctionBegin; 1828 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1829 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1830 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1831 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1832 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1833 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1834 #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1835 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1836 #endif 1837 PetscFunctionReturn(0); 1838 } 1839 1840 #undef __FUNCT__ 1841 #define __FUNCT__ "MatCopy_SeqAIJ" 1842 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1843 { 1844 int ierr; 1845 1846 PetscFunctionBegin; 1847 /* If the two matrices have the same copy implementation, use fast copy. */ 1848 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1849 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1850 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1851 1852 if (a->i[A->m] != b->i[B->m]) { 1853 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1854 } 1855 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m])*sizeof(PetscScalar));CHKERRQ(ierr); 1856 } else { 1857 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1858 } 1859 PetscFunctionReturn(0); 1860 } 1861 1862 #undef __FUNCT__ 1863 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1864 int MatSetUpPreallocation_SeqAIJ(Mat A) 1865 { 1866 int ierr; 1867 1868 PetscFunctionBegin; 1869 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1870 PetscFunctionReturn(0); 1871 } 1872 1873 #undef __FUNCT__ 1874 #define __FUNCT__ "MatGetArray_SeqAIJ" 1875 int MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1876 { 1877 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1878 PetscFunctionBegin; 1879 *array = a->a; 1880 PetscFunctionReturn(0); 1881 } 1882 1883 #undef __FUNCT__ 1884 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1885 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 1886 { 1887 PetscFunctionBegin; 1888 PetscFunctionReturn(0); 1889 } 1890 1891 #undef __FUNCT__ 1892 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1893 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1894 { 1895 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1896 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1897 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1898 PetscScalar *vscale_array; 1899 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1900 Vec w1,w2,w3; 1901 void *fctx = coloring->fctx; 1902 PetscTruth flg; 1903 1904 PetscFunctionBegin; 1905 if (!coloring->w1) { 1906 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1907 PetscLogObjectParent(coloring,coloring->w1); 1908 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1909 PetscLogObjectParent(coloring,coloring->w2); 1910 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1911 PetscLogObjectParent(coloring,coloring->w3); 1912 } 1913 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1914 1915 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1916 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1917 if (flg) { 1918 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1919 } else { 1920 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1921 } 1922 1923 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1924 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1925 1926 /* 1927 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1928 coloring->F for the coarser grids from the finest 1929 */ 1930 if (coloring->F) { 1931 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1932 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1933 if (m1 != m2) { 1934 coloring->F = 0; 1935 } 1936 } 1937 1938 if (coloring->F) { 1939 w1 = coloring->F; 1940 coloring->F = 0; 1941 } else { 1942 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1943 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1944 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 1945 } 1946 1947 /* 1948 Compute all the scale factors and share with other processors 1949 */ 1950 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1951 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1952 for (k=0; k<coloring->ncolors; k++) { 1953 /* 1954 Loop over each column associated with color adding the 1955 perturbation to the vector w3. 1956 */ 1957 for (l=0; l<coloring->ncolumns[k]; l++) { 1958 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1959 dx = xx[col]; 1960 if (dx == 0.0) dx = 1.0; 1961 #if !defined(PETSC_USE_COMPLEX) 1962 if (dx < umin && dx >= 0.0) dx = umin; 1963 else if (dx < 0.0 && dx > -umin) dx = -umin; 1964 #else 1965 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1966 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1967 #endif 1968 dx *= epsilon; 1969 vscale_array[col] = 1.0/dx; 1970 } 1971 } 1972 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1973 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1974 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1975 1976 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1977 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 1978 1979 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 1980 else vscaleforrow = coloring->columnsforrow; 1981 1982 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1983 /* 1984 Loop over each color 1985 */ 1986 for (k=0; k<coloring->ncolors; k++) { 1987 coloring->currentcolor = k; 1988 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 1989 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 1990 /* 1991 Loop over each column associated with color adding the 1992 perturbation to the vector w3. 1993 */ 1994 for (l=0; l<coloring->ncolumns[k]; l++) { 1995 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1996 dx = xx[col]; 1997 if (dx == 0.0) dx = 1.0; 1998 #if !defined(PETSC_USE_COMPLEX) 1999 if (dx < umin && dx >= 0.0) dx = umin; 2000 else if (dx < 0.0 && dx > -umin) dx = -umin; 2001 #else 2002 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2003 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2004 #endif 2005 dx *= epsilon; 2006 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2007 w3_array[col] += dx; 2008 } 2009 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2010 2011 /* 2012 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2013 */ 2014 2015 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2016 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2017 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2018 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2019 2020 /* 2021 Loop over rows of vector, putting results into Jacobian matrix 2022 */ 2023 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2024 for (l=0; l<coloring->nrows[k]; l++) { 2025 row = coloring->rows[k][l]; 2026 col = coloring->columnsforrow[k][l]; 2027 y[row] *= vscale_array[vscaleforrow[k][l]]; 2028 srow = row + start; 2029 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2030 } 2031 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2032 } 2033 coloring->currentcolor = k; 2034 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2035 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2036 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2037 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2038 PetscFunctionReturn(0); 2039 } 2040 2041 #include "petscblaslapack.h" 2042 #undef __FUNCT__ 2043 #define __FUNCT__ "MatAXPY_SeqAIJ" 2044 int MatAXPY_SeqAIJ(const PetscScalar a[],Mat X,Mat Y,MatStructure str) 2045 { 2046 int ierr,one=1,i; 2047 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2048 2049 PetscFunctionBegin; 2050 if (str == SAME_NONZERO_PATTERN) { 2051 BLaxpy_(&x->nz,(PetscScalar*)a,x->a,&one,y->a,&one); 2052 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2053 if (y->xtoy && y->XtoY != X) { 2054 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2055 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2056 } 2057 if (!y->xtoy) { /* get xtoy */ 2058 ierr = MatAXPYGetxtoy_Private(X->m,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2059 y->XtoY = X; 2060 } 2061 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += (*a)*(x->a[i]); 2062 PetscLogInfo(0,"MatAXPY_SeqAIJ: ratio of nnz(X)/nnz(Y): %d/%d = %g\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz); 2063 } else { 2064 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 2065 } 2066 PetscFunctionReturn(0); 2067 } 2068 2069 /* -------------------------------------------------------------------*/ 2070 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2071 MatGetRow_SeqAIJ, 2072 MatRestoreRow_SeqAIJ, 2073 MatMult_SeqAIJ, 2074 /* 4*/ MatMultAdd_SeqAIJ, 2075 MatMultTranspose_SeqAIJ, 2076 MatMultTransposeAdd_SeqAIJ, 2077 MatSolve_SeqAIJ, 2078 MatSolveAdd_SeqAIJ, 2079 MatSolveTranspose_SeqAIJ, 2080 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2081 MatLUFactor_SeqAIJ, 2082 0, 2083 MatRelax_SeqAIJ, 2084 MatTranspose_SeqAIJ, 2085 /*15*/ MatGetInfo_SeqAIJ, 2086 MatEqual_SeqAIJ, 2087 MatGetDiagonal_SeqAIJ, 2088 MatDiagonalScale_SeqAIJ, 2089 MatNorm_SeqAIJ, 2090 /*20*/ 0, 2091 MatAssemblyEnd_SeqAIJ, 2092 MatCompress_SeqAIJ, 2093 MatSetOption_SeqAIJ, 2094 MatZeroEntries_SeqAIJ, 2095 /*25*/ MatZeroRows_SeqAIJ, 2096 MatLUFactorSymbolic_SeqAIJ, 2097 MatLUFactorNumeric_SeqAIJ, 2098 MatCholeskyFactorSymbolic_SeqAIJ, 2099 MatCholeskyFactorNumeric_SeqAIJ, 2100 /*30*/ MatSetUpPreallocation_SeqAIJ, 2101 MatILUFactorSymbolic_SeqAIJ, 2102 MatICCFactorSymbolic_SeqAIJ, 2103 MatGetArray_SeqAIJ, 2104 MatRestoreArray_SeqAIJ, 2105 /*35*/ MatDuplicate_SeqAIJ, 2106 0, 2107 0, 2108 MatILUFactor_SeqAIJ, 2109 0, 2110 /*40*/ MatAXPY_SeqAIJ, 2111 MatGetSubMatrices_SeqAIJ, 2112 MatIncreaseOverlap_SeqAIJ, 2113 MatGetValues_SeqAIJ, 2114 MatCopy_SeqAIJ, 2115 /*45*/ MatPrintHelp_SeqAIJ, 2116 MatScale_SeqAIJ, 2117 0, 2118 0, 2119 MatILUDTFactor_SeqAIJ, 2120 /*50*/ MatGetBlockSize_SeqAIJ, 2121 MatGetRowIJ_SeqAIJ, 2122 MatRestoreRowIJ_SeqAIJ, 2123 MatGetColumnIJ_SeqAIJ, 2124 MatRestoreColumnIJ_SeqAIJ, 2125 /*55*/ MatFDColoringCreate_SeqAIJ, 2126 0, 2127 0, 2128 MatPermute_SeqAIJ, 2129 0, 2130 /*60*/ 0, 2131 MatDestroy_SeqAIJ, 2132 MatView_SeqAIJ, 2133 MatGetPetscMaps_Petsc, 2134 0, 2135 /*65*/ 0, 2136 0, 2137 0, 2138 0, 2139 0, 2140 /*70*/ 0, 2141 0, 2142 MatSetColoring_SeqAIJ, 2143 MatSetValuesAdic_SeqAIJ, 2144 MatSetValuesAdifor_SeqAIJ, 2145 /*75*/ MatFDColoringApply_SeqAIJ, 2146 0, 2147 0, 2148 0, 2149 0, 2150 /*80*/ 0, 2151 0, 2152 0, 2153 0, 2154 /*85*/ MatLoad_SeqAIJ, 2155 MatIsSymmetric_SeqAIJ, 2156 }; 2157 2158 EXTERN_C_BEGIN 2159 #undef __FUNCT__ 2160 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 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 = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2701 "MatSeqAIJSetColumnIndices_SeqAIJ", 2702 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2703 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2704 "MatStoreValues_SeqAIJ", 2705 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2706 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2707 "MatRetrieveValues_SeqAIJ", 2708 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2709 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2710 "MatConvert_SeqAIJ_SeqSBAIJ", 2711 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2712 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2713 "MatConvert_SeqAIJ_SeqBAIJ", 2714 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2715 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2716 "MatIsTranspose_SeqAIJ", 2717 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2718 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2719 "MatSeqAIJSetPreallocation_SeqAIJ", 2720 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2721 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2722 "MatReorderForNonzeroDiagonal_SeqAIJ", 2723 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2724 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2725 "MatAdjustForInodes_SeqAIJ", 2726 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2727 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2728 "MatSeqAIJGetInodeSizes_SeqAIJ", 2729 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2730 #if defined(PETSC_HAVE_MATLAB) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2731 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2732 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2733 #endif 2734 ierr = RegisterApplyPtAPRoutines_Private(B);CHKERRQ(ierr); 2735 PetscFunctionReturn(0); 2736 } 2737 EXTERN_C_END 2738 2739 #undef __FUNCT__ 2740 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2741 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2742 { 2743 Mat C; 2744 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2745 int i,m = A->m,ierr; 2746 size_t len; 2747 2748 PetscFunctionBegin; 2749 *B = 0; 2750 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2751 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 2752 c = (Mat_SeqAIJ*)C->data; 2753 2754 C->factor = A->factor; 2755 c->row = 0; 2756 c->col = 0; 2757 c->icol = 0; 2758 c->keepzeroedrows = a->keepzeroedrows; 2759 C->assembled = PETSC_TRUE; 2760 2761 C->M = A->m; 2762 C->N = A->n; 2763 2764 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2765 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2766 for (i=0; i<m; i++) { 2767 c->imax[i] = a->imax[i]; 2768 c->ilen[i] = a->ilen[i]; 2769 } 2770 2771 /* allocate the matrix space */ 2772 c->singlemalloc = PETSC_TRUE; 2773 len = ((size_t) (m+1))*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2774 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2775 c->j = (int*)(c->a + a->i[m] ); 2776 c->i = c->j + a->i[m]; 2777 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2778 if (m > 0) { 2779 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(int));CHKERRQ(ierr); 2780 if (cpvalues == MAT_COPY_VALUES) { 2781 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2782 } else { 2783 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2784 } 2785 } 2786 2787 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2788 c->sorted = a->sorted; 2789 c->roworiented = a->roworiented; 2790 c->nonew = a->nonew; 2791 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2792 c->saved_values = 0; 2793 c->idiag = 0; 2794 c->ssor = 0; 2795 c->ignorezeroentries = a->ignorezeroentries; 2796 c->freedata = PETSC_TRUE; 2797 2798 if (a->diag) { 2799 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2800 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2801 for (i=0; i<m; i++) { 2802 c->diag[i] = a->diag[i]; 2803 } 2804 } else c->diag = 0; 2805 c->inode.use = a->inode.use; 2806 c->inode.limit = a->inode.limit; 2807 c->inode.max_limit = a->inode.max_limit; 2808 if (a->inode.size){ 2809 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2810 c->inode.node_count = a->inode.node_count; 2811 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2812 } else { 2813 c->inode.size = 0; 2814 c->inode.node_count = 0; 2815 } 2816 c->nz = a->nz; 2817 c->maxnz = a->maxnz; 2818 c->solve_work = 0; 2819 C->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2820 C->preallocated = PETSC_TRUE; 2821 2822 *B = C; 2823 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2824 PetscFunctionReturn(0); 2825 } 2826 2827 #undef __FUNCT__ 2828 #define __FUNCT__ "MatLoad_SeqAIJ" 2829 int MatLoad_SeqAIJ(PetscViewer viewer,const MatType type,Mat *A) 2830 { 2831 Mat_SeqAIJ *a; 2832 Mat B; 2833 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N; 2834 MPI_Comm comm; 2835 2836 PetscFunctionBegin; 2837 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2838 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2839 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2840 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2841 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2842 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2843 M = header[1]; N = header[2]; nz = header[3]; 2844 2845 if (nz < 0) { 2846 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2847 } 2848 2849 /* read in row lengths */ 2850 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2851 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2852 2853 /* create our matrix */ 2854 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2855 ierr = MatSetType(B,type);CHKERRQ(ierr); 2856 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2857 a = (Mat_SeqAIJ*)B->data; 2858 2859 /* read in column indices and adjust for Fortran indexing*/ 2860 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2861 2862 /* read in nonzero values */ 2863 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2864 2865 /* set matrix "i" values */ 2866 a->i[0] = 0; 2867 for (i=1; i<= M; i++) { 2868 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2869 a->ilen[i-1] = rowlengths[i-1]; 2870 } 2871 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2872 2873 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2874 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2875 *A = B; 2876 PetscFunctionReturn(0); 2877 } 2878 2879 #undef __FUNCT__ 2880 #define __FUNCT__ "MatEqual_SeqAIJ" 2881 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2882 { 2883 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2884 int ierr; 2885 2886 PetscFunctionBegin; 2887 /* If the matrix dimensions are not equal,or no of nonzeros */ 2888 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2889 *flg = PETSC_FALSE; 2890 PetscFunctionReturn(0); 2891 } 2892 2893 /* if the a->i are the same */ 2894 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2895 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2896 2897 /* if a->j are the same */ 2898 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2899 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2900 2901 /* if a->a are the same */ 2902 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2903 2904 PetscFunctionReturn(0); 2905 2906 } 2907 2908 #undef __FUNCT__ 2909 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2910 /*@C 2911 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2912 provided by the user. 2913 2914 Coolective on MPI_Comm 2915 2916 Input Parameters: 2917 + comm - must be an MPI communicator of size 1 2918 . m - number of rows 2919 . n - number of columns 2920 . i - row indices 2921 . j - column indices 2922 - a - matrix values 2923 2924 Output Parameter: 2925 . mat - the matrix 2926 2927 Level: intermediate 2928 2929 Notes: 2930 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2931 once the matrix is destroyed 2932 2933 You cannot set new nonzero locations into this matrix, that will generate an error. 2934 2935 The i and j indices are 0 based 2936 2937 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2938 2939 @*/ 2940 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2941 { 2942 int ierr,ii; 2943 Mat_SeqAIJ *aij; 2944 2945 PetscFunctionBegin; 2946 ierr = MatCreate(comm,m,n,m,n,mat);CHKERRQ(ierr); 2947 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 2948 ierr = MatSeqAIJSetPreallocation(*mat,SKIP_ALLOCATION,0);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