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