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