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