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