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