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