1 #define PETSCMAT_DLL 2 3 4 /* 5 Defines the basic matrix operations for the AIJ (compressed row) 6 matrix storage format. 7 */ 8 9 10 #include "../src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 11 #include "petscblaslapack.h" 12 #include "petscbt.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 16 PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 17 { 18 PetscErrorCode ierr; 19 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 20 PetscInt i,*diag, m = Y->rmap->n; 21 MatScalar *aa = aij->a; 22 PetscScalar *v; 23 PetscTruth missing; 24 25 PetscFunctionBegin; 26 if (Y->assembled) { 27 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,PETSC_NULL);CHKERRQ(ierr); 28 if (!missing) { 29 diag = aij->diag; 30 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 31 if (is == INSERT_VALUES) { 32 for (i=0; i<m; i++) { 33 aa[diag[i]] = v[i]; 34 } 35 } else { 36 for (i=0; i<m; i++) { 37 aa[diag[i]] += v[i]; 38 } 39 } 40 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 } 44 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 50 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 51 { 52 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 53 PetscErrorCode ierr; 54 PetscInt i,ishift; 55 56 PetscFunctionBegin; 57 *m = A->rmap->n; 58 if (!ia) PetscFunctionReturn(0); 59 ishift = 0; 60 if (symmetric && !A->structurally_symmetric) { 61 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 62 } else if (oshift == 1) { 63 PetscInt nz = a->i[A->rmap->n]; 64 /* malloc space and add 1 to i and j indices */ 65 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 66 for (i=0; i<A->rmap->n+1; i++) (*ia)[i] = a->i[i] + 1; 67 if (ja) { 68 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 69 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 70 } 71 } else { 72 *ia = a->i; 73 if (ja) *ja = a->j; 74 } 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 80 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 81 { 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 if (!ia) PetscFunctionReturn(0); 86 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 87 ierr = PetscFree(*ia);CHKERRQ(ierr); 88 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 89 } 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 95 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 96 { 97 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 98 PetscErrorCode ierr; 99 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 100 PetscInt nz = a->i[m],row,*jj,mr,col; 101 102 PetscFunctionBegin; 103 *nn = n; 104 if (!ia) PetscFunctionReturn(0); 105 if (symmetric) { 106 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 107 } else { 108 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 109 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 110 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 111 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 112 jj = a->j; 113 for (i=0; i<nz; i++) { 114 collengths[jj[i]]++; 115 } 116 cia[0] = oshift; 117 for (i=0; i<n; i++) { 118 cia[i+1] = cia[i] + collengths[i]; 119 } 120 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 121 jj = a->j; 122 for (row=0; row<m; row++) { 123 mr = a->i[row+1] - a->i[row]; 124 for (i=0; i<mr; i++) { 125 col = *jj++; 126 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 127 } 128 } 129 ierr = PetscFree(collengths);CHKERRQ(ierr); 130 *ia = cia; *ja = cja; 131 } 132 PetscFunctionReturn(0); 133 } 134 135 #undef __FUNCT__ 136 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 137 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 138 { 139 PetscErrorCode ierr; 140 141 PetscFunctionBegin; 142 if (!ia) PetscFunctionReturn(0); 143 144 ierr = PetscFree(*ia);CHKERRQ(ierr); 145 ierr = PetscFree(*ja);CHKERRQ(ierr); 146 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 152 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 153 { 154 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 155 PetscInt *ai = a->i; 156 PetscErrorCode ierr; 157 158 PetscFunctionBegin; 159 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatSetValues_SeqAIJ" 165 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 166 { 167 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 168 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 169 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 170 PetscErrorCode ierr; 171 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 172 MatScalar *ap,value,*aa = a->a; 173 PetscTruth ignorezeroentries = a->ignorezeroentries; 174 PetscTruth roworiented = a->roworiented; 175 176 PetscFunctionBegin; 177 if (v) PetscValidScalarPointer(v,6); 178 for (k=0; k<m; k++) { /* loop over added rows */ 179 row = im[k]; 180 if (row < 0) continue; 181 #if defined(PETSC_USE_DEBUG) 182 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 183 #endif 184 rp = aj + ai[row]; ap = aa + ai[row]; 185 rmax = imax[row]; nrow = ailen[row]; 186 low = 0; 187 high = nrow; 188 for (l=0; l<n; l++) { /* loop over added columns */ 189 if (in[l] < 0) continue; 190 #if defined(PETSC_USE_DEBUG) 191 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 192 #endif 193 col = in[l]; 194 if (v) { 195 if (roworiented) { 196 value = v[l + k*n]; 197 } else { 198 value = v[k + l*m]; 199 } 200 } else { 201 value = 0.; 202 } 203 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 204 205 if (col <= lastcol) low = 0; else high = nrow; 206 lastcol = col; 207 while (high-low > 5) { 208 t = (low+high)/2; 209 if (rp[t] > col) high = t; 210 else low = t; 211 } 212 for (i=low; i<high; i++) { 213 if (rp[i] > col) break; 214 if (rp[i] == col) { 215 if (is == ADD_VALUES) ap[i] += value; 216 else ap[i] = value; 217 low = i + 1; 218 goto noinsert; 219 } 220 } 221 if (value == 0.0 && ignorezeroentries) goto noinsert; 222 if (nonew == 1) goto noinsert; 223 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 224 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 225 N = nrow++ - 1; a->nz++; high++; 226 /* shift up all the later entries in this row */ 227 for (ii=N; ii>=i; ii--) { 228 rp[ii+1] = rp[ii]; 229 ap[ii+1] = ap[ii]; 230 } 231 rp[i] = col; 232 ap[i] = value; 233 low = i + 1; 234 noinsert:; 235 } 236 ailen[row] = nrow; 237 } 238 A->same_nonzero = PETSC_FALSE; 239 PetscFunctionReturn(0); 240 } 241 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "MatGetValues_SeqAIJ" 245 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 246 { 247 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 248 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 249 PetscInt *ai = a->i,*ailen = a->ilen; 250 MatScalar *ap,*aa = a->a; 251 252 PetscFunctionBegin; 253 for (k=0; k<m; k++) { /* loop over rows */ 254 row = im[k]; 255 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 256 if (row >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 257 rp = aj + ai[row]; ap = aa + ai[row]; 258 nrow = ailen[row]; 259 for (l=0; l<n; l++) { /* loop over columns */ 260 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 261 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 262 col = in[l] ; 263 high = nrow; low = 0; /* assume unsorted */ 264 while (high-low > 5) { 265 t = (low+high)/2; 266 if (rp[t] > col) high = t; 267 else low = t; 268 } 269 for (i=low; i<high; i++) { 270 if (rp[i] > col) break; 271 if (rp[i] == col) { 272 *v++ = ap[i]; 273 goto finished; 274 } 275 } 276 *v++ = 0.0; 277 finished:; 278 } 279 } 280 PetscFunctionReturn(0); 281 } 282 283 284 #undef __FUNCT__ 285 #define __FUNCT__ "MatView_SeqAIJ_Binary" 286 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 287 { 288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 289 PetscErrorCode ierr; 290 PetscInt i,*col_lens; 291 int fd; 292 293 PetscFunctionBegin; 294 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 295 ierr = PetscMalloc((4+A->rmap->n)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 296 col_lens[0] = MAT_FILE_CLASSID; 297 col_lens[1] = A->rmap->n; 298 col_lens[2] = A->cmap->n; 299 col_lens[3] = a->nz; 300 301 /* store lengths of each row and write (including header) to file */ 302 for (i=0; i<A->rmap->n; i++) { 303 col_lens[4+i] = a->i[i+1] - a->i[i]; 304 } 305 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 306 ierr = PetscFree(col_lens);CHKERRQ(ierr); 307 308 /* store column indices (zero start index) */ 309 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 310 311 /* store nonzero values */ 312 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 317 318 #undef __FUNCT__ 319 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 320 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 321 { 322 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 323 PetscErrorCode ierr; 324 PetscInt i,j,m = A->rmap->n,shift=0; 325 const char *name; 326 PetscViewerFormat format; 327 328 PetscFunctionBegin; 329 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 330 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 331 if (format == PETSC_VIEWER_ASCII_MATLAB) { 332 PetscInt nofinalvalue = 0; 333 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-!shift)) { 334 nofinalvalue = 1; 335 } 336 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 337 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 338 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 339 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 340 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 341 342 for (i=0; i<m; i++) { 343 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 344 #if defined(PETSC_USE_COMPLEX) 345 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); 346 #else 347 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 348 #endif 349 } 350 } 351 if (nofinalvalue) { 352 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 353 } 354 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 355 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 356 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 357 PetscFunctionReturn(0); 358 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 359 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 360 for (i=0; i<m; i++) { 361 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 362 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 363 #if defined(PETSC_USE_COMPLEX) 364 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 365 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 366 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 367 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 368 } else if (PetscRealPart(a->a[j]) != 0.0) { 369 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 370 } 371 #else 372 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 373 #endif 374 } 375 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 376 } 377 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 378 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 379 PetscInt nzd=0,fshift=1,*sptr; 380 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 381 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 382 for (i=0; i<m; i++) { 383 sptr[i] = nzd+1; 384 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 385 if (a->j[j] >= i) { 386 #if defined(PETSC_USE_COMPLEX) 387 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 388 #else 389 if (a->a[j] != 0.0) nzd++; 390 #endif 391 } 392 } 393 } 394 sptr[m] = nzd+1; 395 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 396 for (i=0; i<m+1; i+=6) { 397 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);} 398 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);} 399 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);} 400 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 401 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 402 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 403 } 404 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 405 ierr = PetscFree(sptr);CHKERRQ(ierr); 406 for (i=0; i<m; i++) { 407 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 408 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 409 } 410 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 411 } 412 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 413 for (i=0; i<m; i++) { 414 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 415 if (a->j[j] >= i) { 416 #if defined(PETSC_USE_COMPLEX) 417 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 418 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 419 } 420 #else 421 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 422 #endif 423 } 424 } 425 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 426 } 427 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 428 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 429 PetscInt cnt = 0,jcnt; 430 PetscScalar value; 431 432 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 433 for (i=0; i<m; i++) { 434 jcnt = 0; 435 for (j=0; j<A->cmap->n; j++) { 436 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 437 value = a->a[cnt++]; 438 jcnt++; 439 } else { 440 value = 0.0; 441 } 442 #if defined(PETSC_USE_COMPLEX) 443 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 444 #else 445 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 446 #endif 447 } 448 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 449 } 450 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 451 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 452 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 453 #if defined(PETSC_USE_COMPLEX) 454 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix complex general\n");CHKERRQ(ierr); 455 #else 456 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix real general\n");CHKERRQ(ierr); 457 #endif 458 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 459 for (i=0; i<m; i++) { 460 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 461 #if defined(PETSC_USE_COMPLEX) 462 if (PetscImaginaryPart(a->a[j]) > 0.0) { 463 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 464 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 465 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G -%G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 466 } else { 467 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 468 } 469 #else 470 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %G\n", i+shift, a->j[j]+shift, a->a[j]);CHKERRQ(ierr); 471 #endif 472 } 473 } 474 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 475 } else { 476 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 477 if (A->factortype){ 478 for (i=0; i<m; i++) { 479 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 480 /* L part */ 481 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 482 #if defined(PETSC_USE_COMPLEX) 483 if (PetscImaginaryPart(a->a[j]) > 0.0) { 484 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 485 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 486 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 487 } else { 488 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 489 } 490 #else 491 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 492 #endif 493 } 494 /* diagonal */ 495 j = a->diag[i]; 496 #if defined(PETSC_USE_COMPLEX) 497 if (PetscImaginaryPart(a->a[j]) > 0.0) { 498 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 499 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 500 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 501 } else { 502 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 503 } 504 #else 505 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 506 #endif 507 508 /* U part */ 509 for (j=a->diag[i+1]+1+shift; j<a->diag[i]+shift; j++) { 510 #if defined(PETSC_USE_COMPLEX) 511 if (PetscImaginaryPart(a->a[j]) > 0.0) { 512 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 513 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 514 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 515 } else { 516 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 517 } 518 #else 519 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 520 #endif 521 } 522 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 523 } 524 } else { 525 for (i=0; i<m; i++) { 526 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 527 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 528 #if defined(PETSC_USE_COMPLEX) 529 if (PetscImaginaryPart(a->a[j]) > 0.0) { 530 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 531 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 532 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 533 } else { 534 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 535 } 536 #else 537 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 538 #endif 539 } 540 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 541 } 542 } 543 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 544 } 545 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 546 PetscFunctionReturn(0); 547 } 548 549 #undef __FUNCT__ 550 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 551 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 552 { 553 Mat A = (Mat) Aa; 554 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 555 PetscErrorCode ierr; 556 PetscInt i,j,m = A->rmap->n,color; 557 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 558 PetscViewer viewer; 559 PetscViewerFormat format; 560 561 PetscFunctionBegin; 562 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 563 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 564 565 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 566 /* loop over matrix elements drawing boxes */ 567 568 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 569 /* Blue for negative, Cyan for zero and Red for positive */ 570 color = PETSC_DRAW_BLUE; 571 for (i=0; i<m; i++) { 572 y_l = m - i - 1.0; y_r = y_l + 1.0; 573 for (j=a->i[i]; j<a->i[i+1]; j++) { 574 x_l = a->j[j] ; x_r = x_l + 1.0; 575 #if defined(PETSC_USE_COMPLEX) 576 if (PetscRealPart(a->a[j]) >= 0.) continue; 577 #else 578 if (a->a[j] >= 0.) continue; 579 #endif 580 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 581 } 582 } 583 color = PETSC_DRAW_CYAN; 584 for (i=0; i<m; i++) { 585 y_l = m - i - 1.0; y_r = y_l + 1.0; 586 for (j=a->i[i]; j<a->i[i+1]; j++) { 587 x_l = a->j[j]; x_r = x_l + 1.0; 588 if (a->a[j] != 0.) continue; 589 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 590 } 591 } 592 color = PETSC_DRAW_RED; 593 for (i=0; i<m; i++) { 594 y_l = m - i - 1.0; y_r = y_l + 1.0; 595 for (j=a->i[i]; j<a->i[i+1]; j++) { 596 x_l = a->j[j]; x_r = x_l + 1.0; 597 #if defined(PETSC_USE_COMPLEX) 598 if (PetscRealPart(a->a[j]) <= 0.) continue; 599 #else 600 if (a->a[j] <= 0.) continue; 601 #endif 602 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 603 } 604 } 605 } else { 606 /* use contour shading to indicate magnitude of values */ 607 /* first determine max of all nonzero values */ 608 PetscInt nz = a->nz,count; 609 PetscDraw popup; 610 PetscReal scale; 611 612 for (i=0; i<nz; i++) { 613 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 614 } 615 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 616 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 617 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 618 count = 0; 619 for (i=0; i<m; i++) { 620 y_l = m - i - 1.0; y_r = y_l + 1.0; 621 for (j=a->i[i]; j<a->i[i+1]; j++) { 622 x_l = a->j[j]; x_r = x_l + 1.0; 623 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 624 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 625 count++; 626 } 627 } 628 } 629 PetscFunctionReturn(0); 630 } 631 632 #undef __FUNCT__ 633 #define __FUNCT__ "MatView_SeqAIJ_Draw" 634 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 635 { 636 PetscErrorCode ierr; 637 PetscDraw draw; 638 PetscReal xr,yr,xl,yl,h,w; 639 PetscTruth isnull; 640 641 PetscFunctionBegin; 642 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 643 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 644 if (isnull) PetscFunctionReturn(0); 645 646 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 647 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 648 xr += w; yr += h; xl = -w; yl = -h; 649 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 650 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 651 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 652 PetscFunctionReturn(0); 653 } 654 655 #undef __FUNCT__ 656 #define __FUNCT__ "MatView_SeqAIJ" 657 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 658 { 659 PetscErrorCode ierr; 660 PetscTruth iascii,isbinary,isdraw; 661 662 PetscFunctionBegin; 663 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 664 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 665 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 666 if (iascii) { 667 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 668 } else if (isbinary) { 669 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 670 } else if (isdraw) { 671 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 672 } else { 673 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 674 } 675 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 676 PetscFunctionReturn(0); 677 } 678 679 #undef __FUNCT__ 680 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 681 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 682 { 683 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 684 PetscErrorCode ierr; 685 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 686 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 687 MatScalar *aa = a->a,*ap; 688 PetscReal ratio=0.6; 689 690 PetscFunctionBegin; 691 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 692 693 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 694 for (i=1; i<m; i++) { 695 /* move each row back by the amount of empty slots (fshift) before it*/ 696 fshift += imax[i-1] - ailen[i-1]; 697 rmax = PetscMax(rmax,ailen[i]); 698 if (fshift) { 699 ip = aj + ai[i] ; 700 ap = aa + ai[i] ; 701 N = ailen[i]; 702 for (j=0; j<N; j++) { 703 ip[j-fshift] = ip[j]; 704 ap[j-fshift] = ap[j]; 705 } 706 } 707 ai[i] = ai[i-1] + ailen[i-1]; 708 } 709 if (m) { 710 fshift += imax[m-1] - ailen[m-1]; 711 ai[m] = ai[m-1] + ailen[m-1]; 712 } 713 /* reset ilen and imax for each row */ 714 for (i=0; i<m; i++) { 715 ailen[i] = imax[i] = ai[i+1] - ai[i]; 716 } 717 a->nz = ai[m]; 718 if (fshift && a->nounused == -1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D, %D unneeded", m, A->cmap->n, fshift); 719 720 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 721 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 722 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 723 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 724 A->info.mallocs += a->reallocs; 725 a->reallocs = 0; 726 A->info.nz_unneeded = (double)fshift; 727 a->rmax = rmax; 728 729 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 730 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 731 A->same_nonzero = PETSC_TRUE; 732 733 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 734 735 a->idiagvalid = PETSC_FALSE; 736 737 #if defined(PETSC_HAVE_CUDA) 738 a->GPUmatrix.resize(m+1,A->cmap->n,a->nz); 739 a->GPUmatrix.row_offsets.assign(a->i,a->i+m+1); 740 a->GPUmatrix.column_indices.assign(a->j,a->j+A->cmap->n); 741 a->GPUmatrix.values.assign(a->a,a->a+a->nz); 742 #endif 743 PetscFunctionReturn(0); 744 } 745 746 #undef __FUNCT__ 747 #define __FUNCT__ "MatRealPart_SeqAIJ" 748 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 749 { 750 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 751 PetscInt i,nz = a->nz; 752 MatScalar *aa = a->a; 753 754 PetscFunctionBegin; 755 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 756 PetscFunctionReturn(0); 757 } 758 759 #undef __FUNCT__ 760 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 761 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 762 { 763 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 764 PetscInt i,nz = a->nz; 765 MatScalar *aa = a->a; 766 767 PetscFunctionBegin; 768 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 769 PetscFunctionReturn(0); 770 } 771 772 #undef __FUNCT__ 773 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 774 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 775 { 776 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 777 PetscErrorCode ierr; 778 779 PetscFunctionBegin; 780 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 781 PetscFunctionReturn(0); 782 } 783 784 #undef __FUNCT__ 785 #define __FUNCT__ "MatDestroy_SeqAIJ" 786 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 787 { 788 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 789 PetscErrorCode ierr; 790 791 PetscFunctionBegin; 792 #if defined(PETSC_USE_LOG) 793 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 794 #endif 795 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 796 if (a->row) { 797 ierr = ISDestroy(a->row);CHKERRQ(ierr); 798 } 799 if (a->col) { 800 ierr = ISDestroy(a->col);CHKERRQ(ierr); 801 } 802 ierr = PetscFree(a->diag);CHKERRQ(ierr); 803 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 804 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 805 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 806 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 807 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 808 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 809 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 810 if (a->XtoY) {ierr = MatDestroy(a->XtoY);CHKERRQ(ierr);} 811 if (a->compressedrow.checked && a->compressedrow.use){ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);} 812 813 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 814 815 ierr = PetscFree(a);CHKERRQ(ierr); 816 817 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 818 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 819 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 820 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 821 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 822 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 823 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqcsrperm_C","",PETSC_NULL);CHKERRQ(ierr); 824 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 825 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 826 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C","",PETSC_NULL);CHKERRQ(ierr); 827 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 828 PetscFunctionReturn(0); 829 } 830 831 #undef __FUNCT__ 832 #define __FUNCT__ "MatSetOption_SeqAIJ" 833 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscTruth flg) 834 { 835 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 836 PetscErrorCode ierr; 837 838 PetscFunctionBegin; 839 switch (op) { 840 case MAT_ROW_ORIENTED: 841 a->roworiented = flg; 842 break; 843 case MAT_KEEP_NONZERO_PATTERN: 844 a->keepnonzeropattern = flg; 845 break; 846 case MAT_NEW_NONZERO_LOCATIONS: 847 a->nonew = (flg ? 0 : 1); 848 break; 849 case MAT_NEW_NONZERO_LOCATION_ERR: 850 a->nonew = (flg ? -1 : 0); 851 break; 852 case MAT_NEW_NONZERO_ALLOCATION_ERR: 853 a->nonew = (flg ? -2 : 0); 854 break; 855 case MAT_UNUSED_NONZERO_LOCATION_ERR: 856 a->nounused = (flg ? -1 : 0); 857 break; 858 case MAT_IGNORE_ZERO_ENTRIES: 859 a->ignorezeroentries = flg; 860 break; 861 case MAT_USE_COMPRESSEDROW: 862 a->compressedrow.use = flg; 863 break; 864 case MAT_SYMMETRIC: 865 case MAT_STRUCTURALLY_SYMMETRIC: 866 case MAT_HERMITIAN: 867 case MAT_SYMMETRY_ETERNAL: 868 case MAT_NEW_DIAGONALS: 869 case MAT_IGNORE_OFF_PROC_ENTRIES: 870 case MAT_USE_HASH_TABLE: 871 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 872 break; 873 case MAT_USE_INODES: 874 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 875 break; 876 default: 877 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 878 } 879 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 880 PetscFunctionReturn(0); 881 } 882 883 #undef __FUNCT__ 884 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 885 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 886 { 887 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 888 PetscErrorCode ierr; 889 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 890 PetscScalar *aa=a->a,*x,zero=0.0; 891 892 PetscFunctionBegin; 893 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 894 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 895 896 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU){ 897 PetscInt *diag=a->diag; 898 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 899 for (i=0; i<n; i++) x[i] = aa[diag[i]]; 900 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 901 PetscFunctionReturn(0); 902 } 903 904 ierr = VecSet(v,zero);CHKERRQ(ierr); 905 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 906 for (i=0; i<n; i++) { 907 nz = ai[i+1] - ai[i]; 908 if (!nz) x[i] = 0.0; 909 for (j=ai[i]; j<ai[i+1]; j++){ 910 if (aj[j] == i) { 911 x[i] = aa[j]; 912 break; 913 } 914 } 915 } 916 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 917 PetscFunctionReturn(0); 918 } 919 920 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 921 #undef __FUNCT__ 922 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 923 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 924 { 925 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 926 PetscScalar *x,*y; 927 PetscErrorCode ierr; 928 PetscInt m = A->rmap->n; 929 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 930 MatScalar *v; 931 PetscScalar alpha; 932 PetscInt n,i,j,*idx,*ii,*ridx=PETSC_NULL; 933 Mat_CompressedRow cprow = a->compressedrow; 934 PetscTruth usecprow = cprow.use; 935 #endif 936 937 PetscFunctionBegin; 938 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 939 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 940 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 941 942 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 943 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 944 #else 945 if (usecprow){ 946 m = cprow.nrows; 947 ii = cprow.i; 948 ridx = cprow.rindex; 949 } else { 950 ii = a->i; 951 } 952 for (i=0; i<m; i++) { 953 idx = a->j + ii[i] ; 954 v = a->a + ii[i] ; 955 n = ii[i+1] - ii[i]; 956 if (usecprow){ 957 alpha = x[ridx[i]]; 958 } else { 959 alpha = x[i]; 960 } 961 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 962 } 963 #endif 964 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 965 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 966 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 967 PetscFunctionReturn(0); 968 } 969 970 #undef __FUNCT__ 971 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 972 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 973 { 974 PetscErrorCode ierr; 975 976 PetscFunctionBegin; 977 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 978 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 979 PetscFunctionReturn(0); 980 } 981 982 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 983 #undef __FUNCT__ 984 #define __FUNCT__ "MatMult_SeqAIJ" 985 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 986 { 987 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 988 PetscScalar *y; 989 const PetscScalar *x; 990 const MatScalar *aa; 991 PetscErrorCode ierr; 992 PetscInt m=A->rmap->n; 993 const PetscInt *aj,*ii,*ridx=PETSC_NULL; 994 PetscInt n,i,nonzerorow=0; 995 PetscScalar sum; 996 PetscTruth usecprow=a->compressedrow.use; 997 998 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 999 #pragma disjoint(*x,*y,*aa) 1000 #endif 1001 1002 PetscFunctionBegin; 1003 #if !defined(PETSC_HAVE_CUDA) 1004 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1005 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1006 #endif 1007 aj = a->j; 1008 aa = a->a; 1009 ii = a->i; 1010 if (usecprow){ /* use compressed row format */ 1011 m = a->compressedrow.nrows; 1012 ii = a->compressedrow.i; 1013 ridx = a->compressedrow.rindex; 1014 for (i=0; i<m; i++){ 1015 n = ii[i+1] - ii[i]; 1016 aj = a->j + ii[i]; 1017 aa = a->a + ii[i]; 1018 sum = 0.0; 1019 nonzerorow += (n>0); 1020 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1021 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1022 y[*ridx++] = sum; 1023 } 1024 } else { /* do not use compressed row format */ 1025 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1026 fortranmultaij_(&m,x,ii,aj,aa,y); 1027 #else 1028 #if defined(PETSC_HAVE_CUDA) 1029 ierr = VecCUDACopyToGPU(xx);CHKERRQ(ierr); 1030 ierr = VecCUDACopyToGPU(yy);CHKERRQ(ierr); 1031 cusp::multiply(a->GPUmatrix,xx->GPUarray,yy->GPUarray); 1032 yy->valid_GPU_array = PETSC_CUDA_GPU; 1033 /* ierr = VecCUDACopyFromGPU(yy);CHKERRQ(ierr); 1034 ierr = PetscPrintf(PETSC_COMM_WORLD, "b's flag is %i\n\n",yy->valid_GPU_array);CHKERRQ(ierr); */ 1035 #else 1036 for (i=0; i<m; i++) { 1037 n = ii[i+1] - ii[i]; 1038 aj = a->j + ii[i]; 1039 aa = a->a + ii[i]; 1040 sum = 0.0; 1041 nonzerorow += (n>0); 1042 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1043 y[i] = sum; 1044 } 1045 #endif 1046 #endif 1047 } 1048 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1049 #if !defined(PETSC_HAVE_CUDA) 1050 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1051 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1052 #endif 1053 PetscFunctionReturn(0); 1054 } 1055 1056 #include "../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h" 1057 #undef __FUNCT__ 1058 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1059 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1060 { 1061 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1062 PetscScalar *x,*y,*z; 1063 const MatScalar *aa; 1064 PetscErrorCode ierr; 1065 PetscInt m = A->rmap->n,*aj,*ii; 1066 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1067 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 1068 PetscScalar sum; 1069 PetscTruth usecprow=a->compressedrow.use; 1070 #endif 1071 1072 PetscFunctionBegin; 1073 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1074 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1075 if (zz != yy) { 1076 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1077 } else { 1078 z = y; 1079 } 1080 1081 aj = a->j; 1082 aa = a->a; 1083 ii = a->i; 1084 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1085 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1086 #else 1087 if (usecprow){ /* use compressed row format */ 1088 if (zz != yy){ 1089 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1090 } 1091 m = a->compressedrow.nrows; 1092 ii = a->compressedrow.i; 1093 ridx = a->compressedrow.rindex; 1094 for (i=0; i<m; i++){ 1095 n = ii[i+1] - ii[i]; 1096 aj = a->j + ii[i]; 1097 aa = a->a + ii[i]; 1098 sum = y[*ridx]; 1099 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 1100 z[*ridx++] = sum; 1101 } 1102 } else { /* do not use compressed row format */ 1103 for (i=0; i<m; i++) { 1104 jrow = ii[i]; 1105 n = ii[i+1] - jrow; 1106 sum = y[i]; 1107 for (j=0; j<n; j++) { 1108 sum += aa[jrow]*x[aj[jrow]]; jrow++; 1109 } 1110 z[i] = sum; 1111 } 1112 } 1113 #endif 1114 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1115 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1116 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1117 if (zz != yy) { 1118 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1119 } 1120 PetscFunctionReturn(0); 1121 } 1122 1123 /* 1124 Adds diagonal pointers to sparse matrix structure. 1125 */ 1126 #undef __FUNCT__ 1127 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1128 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1129 { 1130 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1131 PetscErrorCode ierr; 1132 PetscInt i,j,m = A->rmap->n; 1133 1134 PetscFunctionBegin; 1135 if (!a->diag) { 1136 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1137 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1138 } 1139 for (i=0; i<A->rmap->n; i++) { 1140 a->diag[i] = a->i[i+1]; 1141 for (j=a->i[i]; j<a->i[i+1]; j++) { 1142 if (a->j[j] == i) { 1143 a->diag[i] = j; 1144 break; 1145 } 1146 } 1147 } 1148 PetscFunctionReturn(0); 1149 } 1150 1151 /* 1152 Checks for missing diagonals 1153 */ 1154 #undef __FUNCT__ 1155 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1156 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1157 { 1158 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1159 PetscInt *diag,*jj = a->j,i; 1160 1161 PetscFunctionBegin; 1162 *missing = PETSC_FALSE; 1163 if (A->rmap->n > 0 && !jj) { 1164 *missing = PETSC_TRUE; 1165 if (d) *d = 0; 1166 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1167 } else { 1168 diag = a->diag; 1169 for (i=0; i<A->rmap->n; i++) { 1170 if (jj[diag[i]] != i) { 1171 *missing = PETSC_TRUE; 1172 if (d) *d = i; 1173 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1174 } 1175 } 1176 } 1177 PetscFunctionReturn(0); 1178 } 1179 1180 EXTERN_C_BEGIN 1181 #undef __FUNCT__ 1182 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1183 PetscErrorCode PETSCMAT_DLLEXPORT MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1184 { 1185 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1186 PetscErrorCode ierr; 1187 PetscInt i,*diag,m = A->rmap->n; 1188 MatScalar *v = a->a; 1189 PetscScalar *idiag,*mdiag; 1190 1191 PetscFunctionBegin; 1192 if (a->idiagvalid) PetscFunctionReturn(0); 1193 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1194 diag = a->diag; 1195 if (!a->idiag) { 1196 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1197 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1198 v = a->a; 1199 } 1200 mdiag = a->mdiag; 1201 idiag = a->idiag; 1202 1203 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1204 for (i=0; i<m; i++) { 1205 mdiag[i] = v[diag[i]]; 1206 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1207 idiag[i] = 1.0/v[diag[i]]; 1208 } 1209 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1210 } else { 1211 for (i=0; i<m; i++) { 1212 mdiag[i] = v[diag[i]]; 1213 idiag[i] = omega/(fshift + v[diag[i]]); 1214 } 1215 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1216 } 1217 a->idiagvalid = PETSC_TRUE; 1218 PetscFunctionReturn(0); 1219 } 1220 EXTERN_C_END 1221 1222 #include "../src/mat/impls/aij/seq/ftn-kernels/frelax.h" 1223 #undef __FUNCT__ 1224 #define __FUNCT__ "MatSOR_SeqAIJ" 1225 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1226 { 1227 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1228 PetscScalar *x,d,sum,*t,scale; 1229 const MatScalar *v = a->a,*idiag=0,*mdiag; 1230 const PetscScalar *b, *bs,*xb, *ts; 1231 PetscErrorCode ierr; 1232 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1233 const PetscInt *idx,*diag; 1234 1235 PetscFunctionBegin; 1236 its = its*lits; 1237 1238 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1239 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1240 a->fshift = fshift; 1241 a->omega = omega; 1242 1243 diag = a->diag; 1244 t = a->ssor_work; 1245 idiag = a->idiag; 1246 mdiag = a->mdiag; 1247 1248 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1249 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1250 CHKMEMQ; 1251 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1252 if (flag == SOR_APPLY_UPPER) { 1253 /* apply (U + D/omega) to the vector */ 1254 bs = b; 1255 for (i=0; i<m; i++) { 1256 d = fshift + mdiag[i]; 1257 n = a->i[i+1] - diag[i] - 1; 1258 idx = a->j + diag[i] + 1; 1259 v = a->a + diag[i] + 1; 1260 sum = b[i]*d/omega; 1261 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1262 x[i] = sum; 1263 } 1264 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1265 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1266 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1267 PetscFunctionReturn(0); 1268 } 1269 1270 if (flag == SOR_APPLY_LOWER) { 1271 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1272 } else if (flag & SOR_EISENSTAT) { 1273 /* Let A = L + U + D; where L is lower trianglar, 1274 U is upper triangular, E = D/omega; This routine applies 1275 1276 (L + E)^{-1} A (U + E)^{-1} 1277 1278 to a vector efficiently using Eisenstat's trick. 1279 */ 1280 scale = (2.0/omega) - 1.0; 1281 1282 /* x = (E + U)^{-1} b */ 1283 for (i=m-1; i>=0; i--) { 1284 n = a->i[i+1] - diag[i] - 1; 1285 idx = a->j + diag[i] + 1; 1286 v = a->a + diag[i] + 1; 1287 sum = b[i]; 1288 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1289 x[i] = sum*idiag[i]; 1290 } 1291 1292 /* t = b - (2*E - D)x */ 1293 v = a->a; 1294 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1295 1296 /* t = (E + L)^{-1}t */ 1297 ts = t; 1298 diag = a->diag; 1299 for (i=0; i<m; i++) { 1300 n = diag[i] - a->i[i]; 1301 idx = a->j + a->i[i]; 1302 v = a->a + a->i[i]; 1303 sum = t[i]; 1304 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1305 t[i] = sum*idiag[i]; 1306 /* x = x + t */ 1307 x[i] += t[i]; 1308 } 1309 1310 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1311 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1312 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1313 PetscFunctionReturn(0); 1314 } 1315 if (flag & SOR_ZERO_INITIAL_GUESS) { 1316 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1317 for (i=0; i<m; i++) { 1318 n = diag[i] - a->i[i]; 1319 idx = a->j + a->i[i]; 1320 v = a->a + a->i[i]; 1321 sum = b[i]; 1322 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1323 t[i] = sum; 1324 x[i] = sum*idiag[i]; 1325 } 1326 xb = t; 1327 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1328 } else xb = b; 1329 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1330 for (i=m-1; i>=0; i--) { 1331 n = a->i[i+1] - diag[i] - 1; 1332 idx = a->j + diag[i] + 1; 1333 v = a->a + diag[i] + 1; 1334 sum = xb[i]; 1335 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1336 if (xb == b) { 1337 x[i] = sum*idiag[i]; 1338 } else { 1339 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1340 } 1341 } 1342 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1343 } 1344 its--; 1345 } 1346 while (its--) { 1347 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1348 for (i=0; i<m; i++) { 1349 n = a->i[i+1] - a->i[i]; 1350 idx = a->j + a->i[i]; 1351 v = a->a + a->i[i]; 1352 sum = b[i]; 1353 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1354 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1355 } 1356 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1357 } 1358 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1359 for (i=m-1; i>=0; i--) { 1360 n = a->i[i+1] - a->i[i]; 1361 idx = a->j + a->i[i]; 1362 v = a->a + a->i[i]; 1363 sum = b[i]; 1364 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1365 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1366 } 1367 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1368 } 1369 } 1370 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1371 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1372 CHKMEMQ; PetscFunctionReturn(0); 1373 } 1374 1375 1376 #undef __FUNCT__ 1377 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1378 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1379 { 1380 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1381 1382 PetscFunctionBegin; 1383 info->block_size = 1.0; 1384 info->nz_allocated = (double)a->maxnz; 1385 info->nz_used = (double)a->nz; 1386 info->nz_unneeded = (double)(a->maxnz - a->nz); 1387 info->assemblies = (double)A->num_ass; 1388 info->mallocs = (double)A->info.mallocs; 1389 info->memory = ((PetscObject)A)->mem; 1390 if (A->factortype) { 1391 info->fill_ratio_given = A->info.fill_ratio_given; 1392 info->fill_ratio_needed = A->info.fill_ratio_needed; 1393 info->factor_mallocs = A->info.factor_mallocs; 1394 } else { 1395 info->fill_ratio_given = 0; 1396 info->fill_ratio_needed = 0; 1397 info->factor_mallocs = 0; 1398 } 1399 PetscFunctionReturn(0); 1400 } 1401 1402 #undef __FUNCT__ 1403 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1404 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1405 { 1406 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1407 PetscInt i,m = A->rmap->n - 1,d = 0; 1408 PetscErrorCode ierr; 1409 PetscTruth missing; 1410 1411 PetscFunctionBegin; 1412 if (a->keepnonzeropattern) { 1413 for (i=0; i<N; i++) { 1414 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1415 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1416 } 1417 if (diag != 0.0) { 1418 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1419 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1420 for (i=0; i<N; i++) { 1421 a->a[a->diag[rows[i]]] = diag; 1422 } 1423 } 1424 A->same_nonzero = PETSC_TRUE; 1425 } else { 1426 if (diag != 0.0) { 1427 for (i=0; i<N; i++) { 1428 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1429 if (a->ilen[rows[i]] > 0) { 1430 a->ilen[rows[i]] = 1; 1431 a->a[a->i[rows[i]]] = diag; 1432 a->j[a->i[rows[i]]] = rows[i]; 1433 } else { /* in case row was completely empty */ 1434 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1435 } 1436 } 1437 } else { 1438 for (i=0; i<N; i++) { 1439 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1440 a->ilen[rows[i]] = 0; 1441 } 1442 } 1443 A->same_nonzero = PETSC_FALSE; 1444 } 1445 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1446 PetscFunctionReturn(0); 1447 } 1448 1449 #undef __FUNCT__ 1450 #define __FUNCT__ "MatGetRow_SeqAIJ" 1451 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1452 { 1453 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1454 PetscInt *itmp; 1455 1456 PetscFunctionBegin; 1457 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1458 1459 *nz = a->i[row+1] - a->i[row]; 1460 if (v) *v = a->a + a->i[row]; 1461 if (idx) { 1462 itmp = a->j + a->i[row]; 1463 if (*nz) { 1464 *idx = itmp; 1465 } 1466 else *idx = 0; 1467 } 1468 PetscFunctionReturn(0); 1469 } 1470 1471 /* remove this function? */ 1472 #undef __FUNCT__ 1473 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1474 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1475 { 1476 PetscFunctionBegin; 1477 PetscFunctionReturn(0); 1478 } 1479 1480 #undef __FUNCT__ 1481 #define __FUNCT__ "MatNorm_SeqAIJ" 1482 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1483 { 1484 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1485 MatScalar *v = a->a; 1486 PetscReal sum = 0.0; 1487 PetscErrorCode ierr; 1488 PetscInt i,j; 1489 1490 PetscFunctionBegin; 1491 if (type == NORM_FROBENIUS) { 1492 for (i=0; i<a->nz; i++) { 1493 #if defined(PETSC_USE_COMPLEX) 1494 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1495 #else 1496 sum += (*v)*(*v); v++; 1497 #endif 1498 } 1499 *nrm = sqrt(sum); 1500 } else if (type == NORM_1) { 1501 PetscReal *tmp; 1502 PetscInt *jj = a->j; 1503 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1504 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1505 *nrm = 0.0; 1506 for (j=0; j<a->nz; j++) { 1507 tmp[*jj++] += PetscAbsScalar(*v); v++; 1508 } 1509 for (j=0; j<A->cmap->n; j++) { 1510 if (tmp[j] > *nrm) *nrm = tmp[j]; 1511 } 1512 ierr = PetscFree(tmp);CHKERRQ(ierr); 1513 } else if (type == NORM_INFINITY) { 1514 *nrm = 0.0; 1515 for (j=0; j<A->rmap->n; j++) { 1516 v = a->a + a->i[j]; 1517 sum = 0.0; 1518 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1519 sum += PetscAbsScalar(*v); v++; 1520 } 1521 if (sum > *nrm) *nrm = sum; 1522 } 1523 } else { 1524 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1525 } 1526 PetscFunctionReturn(0); 1527 } 1528 1529 #undef __FUNCT__ 1530 #define __FUNCT__ "MatTranspose_SeqAIJ" 1531 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1532 { 1533 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1534 Mat C; 1535 PetscErrorCode ierr; 1536 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1537 MatScalar *array = a->a; 1538 1539 PetscFunctionBegin; 1540 if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1541 1542 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1543 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1544 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1545 1546 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1547 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1548 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1549 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1550 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1551 ierr = PetscFree(col);CHKERRQ(ierr); 1552 } else { 1553 C = *B; 1554 } 1555 1556 for (i=0; i<m; i++) { 1557 len = ai[i+1]-ai[i]; 1558 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1559 array += len; 1560 aj += len; 1561 } 1562 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1563 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1564 1565 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1566 *B = C; 1567 } else { 1568 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1569 } 1570 PetscFunctionReturn(0); 1571 } 1572 1573 EXTERN_C_BEGIN 1574 #undef __FUNCT__ 1575 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1576 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1577 { 1578 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1579 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1580 MatScalar *va,*vb; 1581 PetscErrorCode ierr; 1582 PetscInt ma,na,mb,nb, i; 1583 1584 PetscFunctionBegin; 1585 bij = (Mat_SeqAIJ *) B->data; 1586 1587 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1588 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1589 if (ma!=nb || na!=mb){ 1590 *f = PETSC_FALSE; 1591 PetscFunctionReturn(0); 1592 } 1593 aii = aij->i; bii = bij->i; 1594 adx = aij->j; bdx = bij->j; 1595 va = aij->a; vb = bij->a; 1596 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1597 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1598 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1599 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1600 1601 *f = PETSC_TRUE; 1602 for (i=0; i<ma; i++) { 1603 while (aptr[i]<aii[i+1]) { 1604 PetscInt idc,idr; 1605 PetscScalar vc,vr; 1606 /* column/row index/value */ 1607 idc = adx[aptr[i]]; 1608 idr = bdx[bptr[idc]]; 1609 vc = va[aptr[i]]; 1610 vr = vb[bptr[idc]]; 1611 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1612 *f = PETSC_FALSE; 1613 goto done; 1614 } else { 1615 aptr[i]++; 1616 if (B || i!=idc) bptr[idc]++; 1617 } 1618 } 1619 } 1620 done: 1621 ierr = PetscFree(aptr);CHKERRQ(ierr); 1622 if (B) { 1623 ierr = PetscFree(bptr);CHKERRQ(ierr); 1624 } 1625 PetscFunctionReturn(0); 1626 } 1627 EXTERN_C_END 1628 1629 EXTERN_C_BEGIN 1630 #undef __FUNCT__ 1631 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1632 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1633 { 1634 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1635 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1636 MatScalar *va,*vb; 1637 PetscErrorCode ierr; 1638 PetscInt ma,na,mb,nb, i; 1639 1640 PetscFunctionBegin; 1641 bij = (Mat_SeqAIJ *) B->data; 1642 1643 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1644 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1645 if (ma!=nb || na!=mb){ 1646 *f = PETSC_FALSE; 1647 PetscFunctionReturn(0); 1648 } 1649 aii = aij->i; bii = bij->i; 1650 adx = aij->j; bdx = bij->j; 1651 va = aij->a; vb = bij->a; 1652 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1653 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1654 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1655 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1656 1657 *f = PETSC_TRUE; 1658 for (i=0; i<ma; i++) { 1659 while (aptr[i]<aii[i+1]) { 1660 PetscInt idc,idr; 1661 PetscScalar vc,vr; 1662 /* column/row index/value */ 1663 idc = adx[aptr[i]]; 1664 idr = bdx[bptr[idc]]; 1665 vc = va[aptr[i]]; 1666 vr = vb[bptr[idc]]; 1667 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1668 *f = PETSC_FALSE; 1669 goto done; 1670 } else { 1671 aptr[i]++; 1672 if (B || i!=idc) bptr[idc]++; 1673 } 1674 } 1675 } 1676 done: 1677 ierr = PetscFree(aptr);CHKERRQ(ierr); 1678 if (B) { 1679 ierr = PetscFree(bptr);CHKERRQ(ierr); 1680 } 1681 PetscFunctionReturn(0); 1682 } 1683 EXTERN_C_END 1684 1685 #undef __FUNCT__ 1686 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1687 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1688 { 1689 PetscErrorCode ierr; 1690 PetscFunctionBegin; 1691 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1692 PetscFunctionReturn(0); 1693 } 1694 1695 #undef __FUNCT__ 1696 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1697 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1698 { 1699 PetscErrorCode ierr; 1700 PetscFunctionBegin; 1701 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1702 PetscFunctionReturn(0); 1703 } 1704 1705 #undef __FUNCT__ 1706 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1707 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1708 { 1709 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1710 PetscScalar *l,*r,x; 1711 MatScalar *v; 1712 PetscErrorCode ierr; 1713 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1714 1715 PetscFunctionBegin; 1716 if (ll) { 1717 /* The local size is used so that VecMPI can be passed to this routine 1718 by MatDiagonalScale_MPIAIJ */ 1719 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1720 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1721 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1722 v = a->a; 1723 for (i=0; i<m; i++) { 1724 x = l[i]; 1725 M = a->i[i+1] - a->i[i]; 1726 for (j=0; j<M; j++) { (*v++) *= x;} 1727 } 1728 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1729 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1730 } 1731 if (rr) { 1732 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1733 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1734 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1735 v = a->a; jj = a->j; 1736 for (i=0; i<nz; i++) { 1737 (*v++) *= r[*jj++]; 1738 } 1739 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1740 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1741 } 1742 PetscFunctionReturn(0); 1743 } 1744 1745 #undef __FUNCT__ 1746 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1747 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1748 { 1749 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1750 PetscErrorCode ierr; 1751 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1752 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1753 const PetscInt *irow,*icol; 1754 PetscInt nrows,ncols; 1755 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1756 MatScalar *a_new,*mat_a; 1757 Mat C; 1758 PetscTruth stride,sorted; 1759 1760 PetscFunctionBegin; 1761 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1762 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1763 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1764 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1765 1766 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1767 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1768 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1769 1770 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1771 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1772 if (stride && step == 1) { 1773 /* special case of contiguous rows */ 1774 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1775 /* loop over new rows determining lens and starting points */ 1776 for (i=0; i<nrows; i++) { 1777 kstart = ai[irow[i]]; 1778 kend = kstart + ailen[irow[i]]; 1779 for (k=kstart; k<kend; k++) { 1780 if (aj[k] >= first) { 1781 starts[i] = k; 1782 break; 1783 } 1784 } 1785 sum = 0; 1786 while (k < kend) { 1787 if (aj[k++] >= first+ncols) break; 1788 sum++; 1789 } 1790 lens[i] = sum; 1791 } 1792 /* create submatrix */ 1793 if (scall == MAT_REUSE_MATRIX) { 1794 PetscInt n_cols,n_rows; 1795 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1796 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1797 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1798 C = *B; 1799 } else { 1800 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1801 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1802 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1803 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1804 } 1805 c = (Mat_SeqAIJ*)C->data; 1806 1807 /* loop over rows inserting into submatrix */ 1808 a_new = c->a; 1809 j_new = c->j; 1810 i_new = c->i; 1811 1812 for (i=0; i<nrows; i++) { 1813 ii = starts[i]; 1814 lensi = lens[i]; 1815 for (k=0; k<lensi; k++) { 1816 *j_new++ = aj[ii+k] - first; 1817 } 1818 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1819 a_new += lensi; 1820 i_new[i+1] = i_new[i] + lensi; 1821 c->ilen[i] = lensi; 1822 } 1823 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1824 } else { 1825 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1826 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1827 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1828 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1829 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1830 /* determine lens of each row */ 1831 for (i=0; i<nrows; i++) { 1832 kstart = ai[irow[i]]; 1833 kend = kstart + a->ilen[irow[i]]; 1834 lens[i] = 0; 1835 for (k=kstart; k<kend; k++) { 1836 if (smap[aj[k]]) { 1837 lens[i]++; 1838 } 1839 } 1840 } 1841 /* Create and fill new matrix */ 1842 if (scall == MAT_REUSE_MATRIX) { 1843 PetscTruth equal; 1844 1845 c = (Mat_SeqAIJ *)((*B)->data); 1846 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1847 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1848 if (!equal) { 1849 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1850 } 1851 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1852 C = *B; 1853 } else { 1854 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1855 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1856 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1857 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1858 } 1859 c = (Mat_SeqAIJ *)(C->data); 1860 for (i=0; i<nrows; i++) { 1861 row = irow[i]; 1862 kstart = ai[row]; 1863 kend = kstart + a->ilen[row]; 1864 mat_i = c->i[i]; 1865 mat_j = c->j + mat_i; 1866 mat_a = c->a + mat_i; 1867 mat_ilen = c->ilen + i; 1868 for (k=kstart; k<kend; k++) { 1869 if ((tcol=smap[a->j[k]])) { 1870 *mat_j++ = tcol - 1; 1871 *mat_a++ = a->a[k]; 1872 (*mat_ilen)++; 1873 1874 } 1875 } 1876 } 1877 /* Free work space */ 1878 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1879 ierr = PetscFree(smap);CHKERRQ(ierr); 1880 ierr = PetscFree(lens);CHKERRQ(ierr); 1881 } 1882 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1883 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1884 1885 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1886 *B = C; 1887 PetscFunctionReturn(0); 1888 } 1889 1890 #undef __FUNCT__ 1891 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1892 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1893 { 1894 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1895 PetscErrorCode ierr; 1896 Mat outA; 1897 PetscTruth row_identity,col_identity; 1898 1899 PetscFunctionBegin; 1900 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1901 1902 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1903 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1904 1905 outA = inA; 1906 outA->factortype = MAT_FACTOR_LU; 1907 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1908 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1909 a->row = row; 1910 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1911 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1912 a->col = col; 1913 1914 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1915 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1916 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1917 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1918 1919 if (!a->solve_work) { /* this matrix may have been factored before */ 1920 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1921 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1922 } 1923 1924 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1925 if (row_identity && col_identity) { 1926 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 1927 } else { 1928 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1929 } 1930 PetscFunctionReturn(0); 1931 } 1932 1933 #undef __FUNCT__ 1934 #define __FUNCT__ "MatScale_SeqAIJ" 1935 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1936 { 1937 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1938 PetscScalar oalpha = alpha; 1939 PetscErrorCode ierr; 1940 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1941 1942 PetscFunctionBegin; 1943 BLASscal_(&bnz,&oalpha,a->a,&one); 1944 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1945 PetscFunctionReturn(0); 1946 } 1947 1948 #undef __FUNCT__ 1949 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1950 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1951 { 1952 PetscErrorCode ierr; 1953 PetscInt i; 1954 1955 PetscFunctionBegin; 1956 if (scall == MAT_INITIAL_MATRIX) { 1957 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1958 } 1959 1960 for (i=0; i<n; i++) { 1961 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1962 } 1963 PetscFunctionReturn(0); 1964 } 1965 1966 #undef __FUNCT__ 1967 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1968 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1969 { 1970 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1971 PetscErrorCode ierr; 1972 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1973 const PetscInt *idx; 1974 PetscInt start,end,*ai,*aj; 1975 PetscBT table; 1976 1977 PetscFunctionBegin; 1978 m = A->rmap->n; 1979 ai = a->i; 1980 aj = a->j; 1981 1982 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1983 1984 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1985 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1986 1987 for (i=0; i<is_max; i++) { 1988 /* Initialize the two local arrays */ 1989 isz = 0; 1990 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1991 1992 /* Extract the indices, assume there can be duplicate entries */ 1993 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1994 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1995 1996 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1997 for (j=0; j<n ; ++j){ 1998 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1999 } 2000 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2001 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 2002 2003 k = 0; 2004 for (j=0; j<ov; j++){ /* for each overlap */ 2005 n = isz; 2006 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 2007 row = nidx[k]; 2008 start = ai[row]; 2009 end = ai[row+1]; 2010 for (l = start; l<end ; l++){ 2011 val = aj[l] ; 2012 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 2013 } 2014 } 2015 } 2016 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 2017 } 2018 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 2019 ierr = PetscFree(nidx);CHKERRQ(ierr); 2020 PetscFunctionReturn(0); 2021 } 2022 2023 /* -------------------------------------------------------------- */ 2024 #undef __FUNCT__ 2025 #define __FUNCT__ "MatPermute_SeqAIJ" 2026 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2027 { 2028 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2029 PetscErrorCode ierr; 2030 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2031 const PetscInt *row,*col; 2032 PetscInt *cnew,j,*lens; 2033 IS icolp,irowp; 2034 PetscInt *cwork = PETSC_NULL; 2035 PetscScalar *vwork = PETSC_NULL; 2036 2037 PetscFunctionBegin; 2038 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2039 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2040 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2041 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2042 2043 /* determine lengths of permuted rows */ 2044 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2045 for (i=0; i<m; i++) { 2046 lens[row[i]] = a->i[i+1] - a->i[i]; 2047 } 2048 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2049 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2050 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2051 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2052 ierr = PetscFree(lens);CHKERRQ(ierr); 2053 2054 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2055 for (i=0; i<m; i++) { 2056 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2057 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2058 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2059 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2060 } 2061 ierr = PetscFree(cnew);CHKERRQ(ierr); 2062 (*B)->assembled = PETSC_FALSE; 2063 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2064 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2065 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2066 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2067 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2068 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2069 PetscFunctionReturn(0); 2070 } 2071 2072 #undef __FUNCT__ 2073 #define __FUNCT__ "MatCopy_SeqAIJ" 2074 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2075 { 2076 PetscErrorCode ierr; 2077 2078 PetscFunctionBegin; 2079 /* If the two matrices have the same copy implementation, use fast copy. */ 2080 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2081 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2082 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2083 2084 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2085 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2086 } 2087 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2088 } else { 2089 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2090 } 2091 PetscFunctionReturn(0); 2092 } 2093 2094 #undef __FUNCT__ 2095 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2096 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2097 { 2098 PetscErrorCode ierr; 2099 2100 PetscFunctionBegin; 2101 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2102 PetscFunctionReturn(0); 2103 } 2104 2105 #undef __FUNCT__ 2106 #define __FUNCT__ "MatGetArray_SeqAIJ" 2107 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2108 { 2109 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2110 PetscFunctionBegin; 2111 *array = a->a; 2112 PetscFunctionReturn(0); 2113 } 2114 2115 #undef __FUNCT__ 2116 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2117 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2118 { 2119 PetscFunctionBegin; 2120 PetscFunctionReturn(0); 2121 } 2122 2123 #undef __FUNCT__ 2124 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2125 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2126 { 2127 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2128 PetscErrorCode ierr; 2129 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2130 PetscScalar dx,*y,*xx,*w3_array; 2131 PetscScalar *vscale_array; 2132 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2133 Vec w1,w2,w3; 2134 void *fctx = coloring->fctx; 2135 PetscTruth flg = PETSC_FALSE; 2136 2137 PetscFunctionBegin; 2138 if (!coloring->w1) { 2139 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2140 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2141 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2142 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2143 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2144 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2145 } 2146 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2147 2148 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2149 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2150 if (flg) { 2151 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2152 } else { 2153 PetscTruth assembled; 2154 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2155 if (assembled) { 2156 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2157 } 2158 } 2159 2160 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2161 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2162 2163 /* 2164 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2165 coloring->F for the coarser grids from the finest 2166 */ 2167 if (coloring->F) { 2168 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2169 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2170 if (m1 != m2) { 2171 coloring->F = 0; 2172 } 2173 } 2174 2175 if (coloring->F) { 2176 w1 = coloring->F; 2177 coloring->F = 0; 2178 } else { 2179 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2180 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2181 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2182 } 2183 2184 /* 2185 Compute all the scale factors and share with other processors 2186 */ 2187 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2188 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2189 for (k=0; k<coloring->ncolors; k++) { 2190 /* 2191 Loop over each column associated with color adding the 2192 perturbation to the vector w3. 2193 */ 2194 for (l=0; l<coloring->ncolumns[k]; l++) { 2195 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2196 dx = xx[col]; 2197 if (dx == 0.0) dx = 1.0; 2198 #if !defined(PETSC_USE_COMPLEX) 2199 if (dx < umin && dx >= 0.0) dx = umin; 2200 else if (dx < 0.0 && dx > -umin) dx = -umin; 2201 #else 2202 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2203 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2204 #endif 2205 dx *= epsilon; 2206 vscale_array[col] = 1.0/dx; 2207 } 2208 } 2209 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2210 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2211 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2212 2213 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2214 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2215 2216 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2217 else vscaleforrow = coloring->columnsforrow; 2218 2219 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2220 /* 2221 Loop over each color 2222 */ 2223 for (k=0; k<coloring->ncolors; k++) { 2224 coloring->currentcolor = k; 2225 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2226 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2227 /* 2228 Loop over each column associated with color adding the 2229 perturbation to the vector w3. 2230 */ 2231 for (l=0; l<coloring->ncolumns[k]; l++) { 2232 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2233 dx = xx[col]; 2234 if (dx == 0.0) dx = 1.0; 2235 #if !defined(PETSC_USE_COMPLEX) 2236 if (dx < umin && dx >= 0.0) dx = umin; 2237 else if (dx < 0.0 && dx > -umin) dx = -umin; 2238 #else 2239 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2240 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2241 #endif 2242 dx *= epsilon; 2243 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2244 w3_array[col] += dx; 2245 } 2246 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2247 2248 /* 2249 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2250 */ 2251 2252 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2253 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2254 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2255 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2256 2257 /* 2258 Loop over rows of vector, putting results into Jacobian matrix 2259 */ 2260 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2261 for (l=0; l<coloring->nrows[k]; l++) { 2262 row = coloring->rows[k][l]; 2263 col = coloring->columnsforrow[k][l]; 2264 y[row] *= vscale_array[vscaleforrow[k][l]]; 2265 srow = row + start; 2266 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2267 } 2268 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2269 } 2270 coloring->currentcolor = k; 2271 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2272 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2273 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2274 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2275 PetscFunctionReturn(0); 2276 } 2277 2278 #undef __FUNCT__ 2279 #define __FUNCT__ "MatAXPY_SeqAIJ" 2280 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2281 { 2282 PetscErrorCode ierr; 2283 PetscInt i; 2284 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2285 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2286 2287 PetscFunctionBegin; 2288 if (str == SAME_NONZERO_PATTERN) { 2289 PetscScalar alpha = a; 2290 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2291 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2292 if (y->xtoy && y->XtoY != X) { 2293 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2294 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2295 } 2296 if (!y->xtoy) { /* get xtoy */ 2297 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2298 y->XtoY = X; 2299 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2300 } 2301 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2302 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2303 } else { 2304 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2305 } 2306 PetscFunctionReturn(0); 2307 } 2308 2309 #undef __FUNCT__ 2310 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2311 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2312 { 2313 PetscErrorCode ierr; 2314 2315 PetscFunctionBegin; 2316 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2317 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2318 PetscFunctionReturn(0); 2319 } 2320 2321 #undef __FUNCT__ 2322 #define __FUNCT__ "MatConjugate_SeqAIJ" 2323 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2324 { 2325 #if defined(PETSC_USE_COMPLEX) 2326 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2327 PetscInt i,nz; 2328 PetscScalar *a; 2329 2330 PetscFunctionBegin; 2331 nz = aij->nz; 2332 a = aij->a; 2333 for (i=0; i<nz; i++) { 2334 a[i] = PetscConj(a[i]); 2335 } 2336 #else 2337 PetscFunctionBegin; 2338 #endif 2339 PetscFunctionReturn(0); 2340 } 2341 2342 #undef __FUNCT__ 2343 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2344 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2345 { 2346 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2347 PetscErrorCode ierr; 2348 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2349 PetscReal atmp; 2350 PetscScalar *x; 2351 MatScalar *aa; 2352 2353 PetscFunctionBegin; 2354 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2355 aa = a->a; 2356 ai = a->i; 2357 aj = a->j; 2358 2359 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2360 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2361 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2362 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2363 for (i=0; i<m; i++) { 2364 ncols = ai[1] - ai[0]; ai++; 2365 x[i] = 0.0; 2366 for (j=0; j<ncols; j++){ 2367 atmp = PetscAbsScalar(*aa); 2368 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2369 aa++; aj++; 2370 } 2371 } 2372 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2373 PetscFunctionReturn(0); 2374 } 2375 2376 #undef __FUNCT__ 2377 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2378 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2379 { 2380 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2381 PetscErrorCode ierr; 2382 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2383 PetscScalar *x; 2384 MatScalar *aa; 2385 2386 PetscFunctionBegin; 2387 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2388 aa = a->a; 2389 ai = a->i; 2390 aj = a->j; 2391 2392 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2393 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2394 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2395 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2396 for (i=0; i<m; i++) { 2397 ncols = ai[1] - ai[0]; ai++; 2398 if (ncols == A->cmap->n) { /* row is dense */ 2399 x[i] = *aa; if (idx) idx[i] = 0; 2400 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2401 x[i] = 0.0; 2402 if (idx) { 2403 idx[i] = 0; /* in case ncols is zero */ 2404 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2405 if (aj[j] > j) { 2406 idx[i] = j; 2407 break; 2408 } 2409 } 2410 } 2411 } 2412 for (j=0; j<ncols; j++){ 2413 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2414 aa++; aj++; 2415 } 2416 } 2417 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2418 PetscFunctionReturn(0); 2419 } 2420 2421 #undef __FUNCT__ 2422 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2423 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2424 { 2425 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2426 PetscErrorCode ierr; 2427 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2428 PetscReal atmp; 2429 PetscScalar *x; 2430 MatScalar *aa; 2431 2432 PetscFunctionBegin; 2433 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2434 aa = a->a; 2435 ai = a->i; 2436 aj = a->j; 2437 2438 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2439 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2440 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2441 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2442 for (i=0; i<m; i++) { 2443 ncols = ai[1] - ai[0]; ai++; 2444 if (ncols) { 2445 /* Get first nonzero */ 2446 for(j = 0; j < ncols; j++) { 2447 atmp = PetscAbsScalar(aa[j]); 2448 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2449 } 2450 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2451 } else { 2452 x[i] = 0.0; if (idx) idx[i] = 0; 2453 } 2454 for(j = 0; j < ncols; j++) { 2455 atmp = PetscAbsScalar(*aa); 2456 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2457 aa++; aj++; 2458 } 2459 } 2460 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2461 PetscFunctionReturn(0); 2462 } 2463 2464 #undef __FUNCT__ 2465 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2466 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2467 { 2468 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2469 PetscErrorCode ierr; 2470 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2471 PetscScalar *x; 2472 MatScalar *aa; 2473 2474 PetscFunctionBegin; 2475 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2476 aa = a->a; 2477 ai = a->i; 2478 aj = a->j; 2479 2480 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2481 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2482 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2483 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2484 for (i=0; i<m; i++) { 2485 ncols = ai[1] - ai[0]; ai++; 2486 if (ncols == A->cmap->n) { /* row is dense */ 2487 x[i] = *aa; if (idx) idx[i] = 0; 2488 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2489 x[i] = 0.0; 2490 if (idx) { /* find first implicit 0.0 in the row */ 2491 idx[i] = 0; /* in case ncols is zero */ 2492 for (j=0;j<ncols;j++) { 2493 if (aj[j] > j) { 2494 idx[i] = j; 2495 break; 2496 } 2497 } 2498 } 2499 } 2500 for (j=0; j<ncols; j++){ 2501 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2502 aa++; aj++; 2503 } 2504 } 2505 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2506 PetscFunctionReturn(0); 2507 } 2508 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2509 /* -------------------------------------------------------------------*/ 2510 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2511 MatGetRow_SeqAIJ, 2512 MatRestoreRow_SeqAIJ, 2513 MatMult_SeqAIJ, 2514 /* 4*/ MatMultAdd_SeqAIJ, 2515 MatMultTranspose_SeqAIJ, 2516 MatMultTransposeAdd_SeqAIJ, 2517 0, 2518 0, 2519 0, 2520 /*10*/ 0, 2521 MatLUFactor_SeqAIJ, 2522 0, 2523 MatSOR_SeqAIJ, 2524 MatTranspose_SeqAIJ, 2525 /*15*/ MatGetInfo_SeqAIJ, 2526 MatEqual_SeqAIJ, 2527 MatGetDiagonal_SeqAIJ, 2528 MatDiagonalScale_SeqAIJ, 2529 MatNorm_SeqAIJ, 2530 /*20*/ 0, 2531 MatAssemblyEnd_SeqAIJ, 2532 MatSetOption_SeqAIJ, 2533 MatZeroEntries_SeqAIJ, 2534 /*24*/ MatZeroRows_SeqAIJ, 2535 0, 2536 0, 2537 0, 2538 0, 2539 /*29*/ MatSetUpPreallocation_SeqAIJ, 2540 0, 2541 0, 2542 MatGetArray_SeqAIJ, 2543 MatRestoreArray_SeqAIJ, 2544 /*34*/ MatDuplicate_SeqAIJ, 2545 0, 2546 0, 2547 MatILUFactor_SeqAIJ, 2548 0, 2549 /*39*/ MatAXPY_SeqAIJ, 2550 MatGetSubMatrices_SeqAIJ, 2551 MatIncreaseOverlap_SeqAIJ, 2552 MatGetValues_SeqAIJ, 2553 MatCopy_SeqAIJ, 2554 /*44*/ MatGetRowMax_SeqAIJ, 2555 MatScale_SeqAIJ, 2556 0, 2557 MatDiagonalSet_SeqAIJ, 2558 0, 2559 /*49*/ MatSetBlockSize_SeqAIJ, 2560 MatGetRowIJ_SeqAIJ, 2561 MatRestoreRowIJ_SeqAIJ, 2562 MatGetColumnIJ_SeqAIJ, 2563 MatRestoreColumnIJ_SeqAIJ, 2564 /*54*/ MatFDColoringCreate_SeqAIJ, 2565 0, 2566 0, 2567 MatPermute_SeqAIJ, 2568 0, 2569 /*59*/ 0, 2570 MatDestroy_SeqAIJ, 2571 MatView_SeqAIJ, 2572 0, 2573 0, 2574 /*64*/ 0, 2575 0, 2576 0, 2577 0, 2578 0, 2579 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2580 MatGetRowMinAbs_SeqAIJ, 2581 0, 2582 MatSetColoring_SeqAIJ, 2583 #if defined(PETSC_HAVE_ADIC) 2584 MatSetValuesAdic_SeqAIJ, 2585 #else 2586 0, 2587 #endif 2588 /*74*/ MatSetValuesAdifor_SeqAIJ, 2589 MatFDColoringApply_AIJ, 2590 0, 2591 0, 2592 0, 2593 /*79*/ 0, 2594 0, 2595 0, 2596 0, 2597 MatLoad_SeqAIJ, 2598 /*84*/ MatIsSymmetric_SeqAIJ, 2599 MatIsHermitian_SeqAIJ, 2600 0, 2601 0, 2602 0, 2603 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2604 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2605 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2606 MatPtAP_Basic, 2607 MatPtAPSymbolic_SeqAIJ, 2608 /*94*/ MatPtAPNumeric_SeqAIJ, 2609 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2610 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2611 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2612 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2613 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2614 0, 2615 0, 2616 MatConjugate_SeqAIJ, 2617 0, 2618 /*104*/MatSetValuesRow_SeqAIJ, 2619 MatRealPart_SeqAIJ, 2620 MatImaginaryPart_SeqAIJ, 2621 0, 2622 0, 2623 /*109*/0, 2624 0, 2625 MatGetRowMin_SeqAIJ, 2626 0, 2627 MatMissingDiagonal_SeqAIJ, 2628 /*114*/0, 2629 0, 2630 0, 2631 0, 2632 0, 2633 /*119*/0, 2634 0, 2635 0, 2636 0, 2637 MatLoadnew_SeqAIJ 2638 }; 2639 2640 EXTERN_C_BEGIN 2641 #undef __FUNCT__ 2642 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2643 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2644 { 2645 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2646 PetscInt i,nz,n; 2647 2648 PetscFunctionBegin; 2649 2650 nz = aij->maxnz; 2651 n = mat->rmap->n; 2652 for (i=0; i<nz; i++) { 2653 aij->j[i] = indices[i]; 2654 } 2655 aij->nz = nz; 2656 for (i=0; i<n; i++) { 2657 aij->ilen[i] = aij->imax[i]; 2658 } 2659 2660 PetscFunctionReturn(0); 2661 } 2662 EXTERN_C_END 2663 2664 #undef __FUNCT__ 2665 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2666 /*@ 2667 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2668 in the matrix. 2669 2670 Input Parameters: 2671 + mat - the SeqAIJ matrix 2672 - indices - the column indices 2673 2674 Level: advanced 2675 2676 Notes: 2677 This can be called if you have precomputed the nonzero structure of the 2678 matrix and want to provide it to the matrix object to improve the performance 2679 of the MatSetValues() operation. 2680 2681 You MUST have set the correct numbers of nonzeros per row in the call to 2682 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2683 2684 MUST be called before any calls to MatSetValues(); 2685 2686 The indices should start with zero, not one. 2687 2688 @*/ 2689 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2690 { 2691 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2692 2693 PetscFunctionBegin; 2694 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2695 PetscValidPointer(indices,2); 2696 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2697 if (f) { 2698 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2699 } else { 2700 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2701 } 2702 PetscFunctionReturn(0); 2703 } 2704 2705 /* ----------------------------------------------------------------------------------------*/ 2706 2707 EXTERN_C_BEGIN 2708 #undef __FUNCT__ 2709 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2710 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2711 { 2712 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2713 PetscErrorCode ierr; 2714 size_t nz = aij->i[mat->rmap->n]; 2715 2716 PetscFunctionBegin; 2717 if (aij->nonew != 1) { 2718 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2719 } 2720 2721 /* allocate space for values if not already there */ 2722 if (!aij->saved_values) { 2723 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2724 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2725 } 2726 2727 /* copy values over */ 2728 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2729 PetscFunctionReturn(0); 2730 } 2731 EXTERN_C_END 2732 2733 #undef __FUNCT__ 2734 #define __FUNCT__ "MatStoreValues" 2735 /*@ 2736 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2737 example, reuse of the linear part of a Jacobian, while recomputing the 2738 nonlinear portion. 2739 2740 Collect on Mat 2741 2742 Input Parameters: 2743 . mat - the matrix (currently only AIJ matrices support this option) 2744 2745 Level: advanced 2746 2747 Common Usage, with SNESSolve(): 2748 $ Create Jacobian matrix 2749 $ Set linear terms into matrix 2750 $ Apply boundary conditions to matrix, at this time matrix must have 2751 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2752 $ boundary conditions again will not change the nonzero structure 2753 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2754 $ ierr = MatStoreValues(mat); 2755 $ Call SNESSetJacobian() with matrix 2756 $ In your Jacobian routine 2757 $ ierr = MatRetrieveValues(mat); 2758 $ Set nonlinear terms in matrix 2759 2760 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2761 $ // build linear portion of Jacobian 2762 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2763 $ ierr = MatStoreValues(mat); 2764 $ loop over nonlinear iterations 2765 $ ierr = MatRetrieveValues(mat); 2766 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2767 $ // call MatAssemblyBegin/End() on matrix 2768 $ Solve linear system with Jacobian 2769 $ endloop 2770 2771 Notes: 2772 Matrix must already be assemblied before calling this routine 2773 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2774 calling this routine. 2775 2776 When this is called multiple times it overwrites the previous set of stored values 2777 and does not allocated additional space. 2778 2779 .seealso: MatRetrieveValues() 2780 2781 @*/ 2782 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2783 { 2784 PetscErrorCode ierr,(*f)(Mat); 2785 2786 PetscFunctionBegin; 2787 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2788 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2789 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2790 2791 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2792 if (f) { 2793 ierr = (*f)(mat);CHKERRQ(ierr); 2794 } else { 2795 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2796 } 2797 PetscFunctionReturn(0); 2798 } 2799 2800 EXTERN_C_BEGIN 2801 #undef __FUNCT__ 2802 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2803 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2804 { 2805 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2806 PetscErrorCode ierr; 2807 PetscInt nz = aij->i[mat->rmap->n]; 2808 2809 PetscFunctionBegin; 2810 if (aij->nonew != 1) { 2811 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2812 } 2813 if (!aij->saved_values) { 2814 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2815 } 2816 /* copy values over */ 2817 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2818 PetscFunctionReturn(0); 2819 } 2820 EXTERN_C_END 2821 2822 #undef __FUNCT__ 2823 #define __FUNCT__ "MatRetrieveValues" 2824 /*@ 2825 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2826 example, reuse of the linear part of a Jacobian, while recomputing the 2827 nonlinear portion. 2828 2829 Collect on Mat 2830 2831 Input Parameters: 2832 . mat - the matrix (currently on AIJ matrices support this option) 2833 2834 Level: advanced 2835 2836 .seealso: MatStoreValues() 2837 2838 @*/ 2839 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2840 { 2841 PetscErrorCode ierr,(*f)(Mat); 2842 2843 PetscFunctionBegin; 2844 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2845 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2846 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2847 2848 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2849 if (f) { 2850 ierr = (*f)(mat);CHKERRQ(ierr); 2851 } else { 2852 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2853 } 2854 PetscFunctionReturn(0); 2855 } 2856 2857 2858 /* --------------------------------------------------------------------------------*/ 2859 #undef __FUNCT__ 2860 #define __FUNCT__ "MatCreateSeqAIJ" 2861 /*@C 2862 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2863 (the default parallel PETSc format). For good matrix assembly performance 2864 the user should preallocate the matrix storage by setting the parameter nz 2865 (or the array nnz). By setting these parameters accurately, performance 2866 during matrix assembly can be increased by more than a factor of 50. 2867 2868 Collective on MPI_Comm 2869 2870 Input Parameters: 2871 + comm - MPI communicator, set to PETSC_COMM_SELF 2872 . m - number of rows 2873 . n - number of columns 2874 . nz - number of nonzeros per row (same for all rows) 2875 - nnz - array containing the number of nonzeros in the various rows 2876 (possibly different for each row) or PETSC_NULL 2877 2878 Output Parameter: 2879 . A - the matrix 2880 2881 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2882 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2883 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2884 2885 Notes: 2886 If nnz is given then nz is ignored 2887 2888 The AIJ format (also called the Yale sparse matrix format or 2889 compressed row storage), is fully compatible with standard Fortran 77 2890 storage. That is, the stored row and column indices can begin at 2891 either one (as in Fortran) or zero. See the users' manual for details. 2892 2893 Specify the preallocated storage with either nz or nnz (not both). 2894 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2895 allocation. For large problems you MUST preallocate memory or you 2896 will get TERRIBLE performance, see the users' manual chapter on matrices. 2897 2898 By default, this format uses inodes (identical nodes) when possible, to 2899 improve numerical efficiency of matrix-vector products and solves. We 2900 search for consecutive rows with the same nonzero structure, thereby 2901 reusing matrix information to achieve increased efficiency. 2902 2903 Options Database Keys: 2904 + -mat_no_inode - Do not use inodes 2905 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2906 - -mat_aij_oneindex - Internally use indexing starting at 1 2907 rather than 0. Note that when calling MatSetValues(), 2908 the user still MUST index entries starting at 0! 2909 2910 Level: intermediate 2911 2912 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2913 2914 @*/ 2915 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2916 { 2917 PetscErrorCode ierr; 2918 2919 PetscFunctionBegin; 2920 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2921 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2922 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2923 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2924 PetscFunctionReturn(0); 2925 } 2926 2927 #undef __FUNCT__ 2928 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2929 /*@C 2930 MatSeqAIJSetPreallocation - For good matrix assembly performance 2931 the user should preallocate the matrix storage by setting the parameter nz 2932 (or the array nnz). By setting these parameters accurately, performance 2933 during matrix assembly can be increased by more than a factor of 50. 2934 2935 Collective on MPI_Comm 2936 2937 Input Parameters: 2938 + B - The matrix-free 2939 . nz - number of nonzeros per row (same for all rows) 2940 - nnz - array containing the number of nonzeros in the various rows 2941 (possibly different for each row) or PETSC_NULL 2942 2943 Notes: 2944 If nnz is given then nz is ignored 2945 2946 The AIJ format (also called the Yale sparse matrix format or 2947 compressed row storage), is fully compatible with standard Fortran 77 2948 storage. That is, the stored row and column indices can begin at 2949 either one (as in Fortran) or zero. See the users' manual for details. 2950 2951 Specify the preallocated storage with either nz or nnz (not both). 2952 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2953 allocation. For large problems you MUST preallocate memory or you 2954 will get TERRIBLE performance, see the users' manual chapter on matrices. 2955 2956 You can call MatGetInfo() to get information on how effective the preallocation was; 2957 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2958 You can also run with the option -info and look for messages with the string 2959 malloc in them to see if additional memory allocation was needed. 2960 2961 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2962 entries or columns indices 2963 2964 By default, this format uses inodes (identical nodes) when possible, to 2965 improve numerical efficiency of matrix-vector products and solves. We 2966 search for consecutive rows with the same nonzero structure, thereby 2967 reusing matrix information to achieve increased efficiency. 2968 2969 Options Database Keys: 2970 + -mat_no_inode - Do not use inodes 2971 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2972 - -mat_aij_oneindex - Internally use indexing starting at 1 2973 rather than 0. Note that when calling MatSetValues(), 2974 the user still MUST index entries starting at 0! 2975 2976 Level: intermediate 2977 2978 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2979 2980 @*/ 2981 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2982 { 2983 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2984 2985 PetscFunctionBegin; 2986 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2987 if (f) { 2988 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2989 } 2990 PetscFunctionReturn(0); 2991 } 2992 2993 EXTERN_C_BEGIN 2994 #undef __FUNCT__ 2995 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2996 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2997 { 2998 Mat_SeqAIJ *b; 2999 PetscTruth skipallocation = PETSC_FALSE; 3000 PetscErrorCode ierr; 3001 PetscInt i; 3002 3003 PetscFunctionBegin; 3004 3005 if (nz == MAT_SKIP_ALLOCATION) { 3006 skipallocation = PETSC_TRUE; 3007 nz = 0; 3008 } 3009 3010 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3011 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3012 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3013 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3014 3015 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3016 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3017 if (nnz) { 3018 for (i=0; i<B->rmap->n; i++) { 3019 if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 3020 if (nnz[i] > B->cmap->n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->cmap->n); 3021 } 3022 } 3023 3024 B->preallocated = PETSC_TRUE; 3025 b = (Mat_SeqAIJ*)B->data; 3026 3027 if (!skipallocation) { 3028 if (!b->imax) { 3029 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3030 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3031 } 3032 if (!nnz) { 3033 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3034 else if (nz <= 0) nz = 1; 3035 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3036 nz = nz*B->rmap->n; 3037 } else { 3038 nz = 0; 3039 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3040 } 3041 /* b->ilen will count nonzeros in each row so far. */ 3042 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3043 3044 /* allocate the matrix space */ 3045 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3046 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3047 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3048 b->i[0] = 0; 3049 for (i=1; i<B->rmap->n+1; i++) { 3050 b->i[i] = b->i[i-1] + b->imax[i-1]; 3051 } 3052 b->singlemalloc = PETSC_TRUE; 3053 b->free_a = PETSC_TRUE; 3054 b->free_ij = PETSC_TRUE; 3055 } else { 3056 b->free_a = PETSC_FALSE; 3057 b->free_ij = PETSC_FALSE; 3058 } 3059 3060 b->nz = 0; 3061 b->maxnz = nz; 3062 B->info.nz_unneeded = (double)b->maxnz; 3063 PetscFunctionReturn(0); 3064 } 3065 EXTERN_C_END 3066 3067 #undef __FUNCT__ 3068 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3069 /*@ 3070 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3071 3072 Input Parameters: 3073 + B - the matrix 3074 . i - the indices into j for the start of each row (starts with zero) 3075 . j - the column indices for each row (starts with zero) these must be sorted for each row 3076 - v - optional values in the matrix 3077 3078 Level: developer 3079 3080 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3081 3082 .keywords: matrix, aij, compressed row, sparse, sequential 3083 3084 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3085 @*/ 3086 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3087 { 3088 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3089 PetscErrorCode ierr; 3090 3091 PetscFunctionBegin; 3092 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3093 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3094 if (f) { 3095 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3096 } 3097 PetscFunctionReturn(0); 3098 } 3099 3100 EXTERN_C_BEGIN 3101 #undef __FUNCT__ 3102 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3103 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3104 { 3105 PetscInt i; 3106 PetscInt m,n; 3107 PetscInt nz; 3108 PetscInt *nnz, nz_max = 0; 3109 PetscScalar *values; 3110 PetscErrorCode ierr; 3111 3112 PetscFunctionBegin; 3113 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3114 3115 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3116 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3117 for(i = 0; i < m; i++) { 3118 nz = Ii[i+1]- Ii[i]; 3119 nz_max = PetscMax(nz_max, nz); 3120 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3121 nnz[i] = nz; 3122 } 3123 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3124 ierr = PetscFree(nnz);CHKERRQ(ierr); 3125 3126 if (v) { 3127 values = (PetscScalar*) v; 3128 } else { 3129 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3130 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3131 } 3132 3133 for(i = 0; i < m; i++) { 3134 nz = Ii[i+1] - Ii[i]; 3135 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3136 } 3137 3138 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3139 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3140 3141 if (!v) { 3142 ierr = PetscFree(values);CHKERRQ(ierr); 3143 } 3144 PetscFunctionReturn(0); 3145 } 3146 EXTERN_C_END 3147 3148 #include "../src/mat/impls/dense/seq/dense.h" 3149 #include "private/petscaxpy.h" 3150 3151 #undef __FUNCT__ 3152 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3153 /* 3154 Computes (B'*A')' since computing B*A directly is untenable 3155 3156 n p p 3157 ( ) ( ) ( ) 3158 m ( A ) * n ( B ) = m ( C ) 3159 ( ) ( ) ( ) 3160 3161 */ 3162 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3163 { 3164 PetscErrorCode ierr; 3165 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3166 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3167 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3168 PetscInt i,n,m,q,p; 3169 const PetscInt *ii,*idx; 3170 const PetscScalar *b,*a,*a_q; 3171 PetscScalar *c,*c_q; 3172 3173 PetscFunctionBegin; 3174 m = A->rmap->n; 3175 n = A->cmap->n; 3176 p = B->cmap->n; 3177 a = sub_a->v; 3178 b = sub_b->a; 3179 c = sub_c->v; 3180 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3181 3182 ii = sub_b->i; 3183 idx = sub_b->j; 3184 for (i=0; i<n; i++) { 3185 q = ii[i+1] - ii[i]; 3186 while (q-->0) { 3187 c_q = c + m*(*idx); 3188 a_q = a + m*i; 3189 PetscAXPY(c_q,*b,a_q,m); 3190 idx++; 3191 b++; 3192 } 3193 } 3194 PetscFunctionReturn(0); 3195 } 3196 3197 #undef __FUNCT__ 3198 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3199 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3200 { 3201 PetscErrorCode ierr; 3202 PetscInt m=A->rmap->n,n=B->cmap->n; 3203 Mat Cmat; 3204 3205 PetscFunctionBegin; 3206 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 3207 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3208 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3209 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3210 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3211 Cmat->assembled = PETSC_TRUE; 3212 *C = Cmat; 3213 PetscFunctionReturn(0); 3214 } 3215 3216 /* ----------------------------------------------------------------*/ 3217 #undef __FUNCT__ 3218 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3219 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3220 { 3221 PetscErrorCode ierr; 3222 3223 PetscFunctionBegin; 3224 if (scall == MAT_INITIAL_MATRIX){ 3225 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3226 } 3227 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3228 PetscFunctionReturn(0); 3229 } 3230 3231 3232 /*MC 3233 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3234 based on compressed sparse row format. 3235 3236 Options Database Keys: 3237 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3238 3239 Level: beginner 3240 3241 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3242 M*/ 3243 3244 EXTERN_C_BEGIN 3245 #if defined(PETSC_HAVE_PASTIX) 3246 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3247 #endif 3248 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3249 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3250 #endif 3251 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3252 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3253 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3254 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3255 #if defined(PETSC_HAVE_MUMPS) 3256 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3257 #endif 3258 #if defined(PETSC_HAVE_SUPERLU) 3259 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3260 #endif 3261 #if defined(PETSC_HAVE_SUPERLU_DIST) 3262 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3263 #endif 3264 #if defined(PETSC_HAVE_SPOOLES) 3265 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3266 #endif 3267 #if defined(PETSC_HAVE_UMFPACK) 3268 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3269 #endif 3270 #if defined(PETSC_HAVE_CHOLMOD) 3271 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3272 #endif 3273 #if defined(PETSC_HAVE_LUSOL) 3274 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3275 #endif 3276 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3277 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3278 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3279 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3280 #endif 3281 EXTERN_C_END 3282 3283 3284 EXTERN_C_BEGIN 3285 #undef __FUNCT__ 3286 #define __FUNCT__ "MatCreate_SeqAIJ" 3287 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3288 { 3289 Mat_SeqAIJ *b; 3290 PetscErrorCode ierr; 3291 PetscMPIInt size; 3292 3293 PetscFunctionBegin; 3294 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3295 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3296 3297 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3298 B->data = (void*)b; 3299 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3300 B->mapping = 0; 3301 b->row = 0; 3302 b->col = 0; 3303 b->icol = 0; 3304 b->reallocs = 0; 3305 b->ignorezeroentries = PETSC_FALSE; 3306 b->roworiented = PETSC_TRUE; 3307 b->nonew = 0; 3308 b->diag = 0; 3309 b->solve_work = 0; 3310 B->spptr = 0; 3311 b->saved_values = 0; 3312 b->idiag = 0; 3313 b->mdiag = 0; 3314 b->ssor_work = 0; 3315 b->omega = 1.0; 3316 b->fshift = 0.0; 3317 b->idiagvalid = PETSC_FALSE; 3318 b->keepnonzeropattern = PETSC_FALSE; 3319 b->xtoy = 0; 3320 b->XtoY = 0; 3321 b->compressedrow.use = PETSC_FALSE; 3322 b->compressedrow.nrows = B->rmap->n; 3323 b->compressedrow.i = PETSC_NULL; 3324 b->compressedrow.rindex = PETSC_NULL; 3325 b->compressedrow.checked = PETSC_FALSE; 3326 B->same_nonzero = PETSC_FALSE; 3327 3328 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3329 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3330 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3331 "MatGetFactor_seqaij_matlab", 3332 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3333 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3334 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3335 #endif 3336 #if defined(PETSC_HAVE_PASTIX) 3337 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3338 "MatGetFactor_seqaij_pastix", 3339 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3340 #endif 3341 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3342 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3343 "MatGetFactor_seqaij_essl", 3344 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3345 #endif 3346 #if defined(PETSC_HAVE_SUPERLU) 3347 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3348 "MatGetFactor_seqaij_superlu", 3349 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3350 #endif 3351 #if defined(PETSC_HAVE_SUPERLU_DIST) 3352 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3353 "MatGetFactor_seqaij_superlu_dist", 3354 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3355 #endif 3356 #if defined(PETSC_HAVE_SPOOLES) 3357 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3358 "MatGetFactor_seqaij_spooles", 3359 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3360 #endif 3361 #if defined(PETSC_HAVE_MUMPS) 3362 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3363 "MatGetFactor_aij_mumps", 3364 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3365 #endif 3366 #if defined(PETSC_HAVE_UMFPACK) 3367 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3368 "MatGetFactor_seqaij_umfpack", 3369 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3370 #endif 3371 #if defined(PETSC_HAVE_CHOLMOD) 3372 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3373 "MatGetFactor_seqaij_cholmod", 3374 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3375 #endif 3376 #if defined(PETSC_HAVE_LUSOL) 3377 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3378 "MatGetFactor_seqaij_lusol", 3379 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3380 #endif 3381 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3382 "MatGetFactor_seqaij_petsc", 3383 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3384 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3385 "MatGetFactorAvailable_seqaij_petsc", 3386 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3387 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3388 "MatGetFactor_seqaij_bas", 3389 MatGetFactor_seqaij_bas); 3390 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3391 "MatSeqAIJSetColumnIndices_SeqAIJ", 3392 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3393 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3394 "MatStoreValues_SeqAIJ", 3395 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3396 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3397 "MatRetrieveValues_SeqAIJ", 3398 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3399 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3400 "MatConvert_SeqAIJ_SeqSBAIJ", 3401 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3402 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3403 "MatConvert_SeqAIJ_SeqBAIJ", 3404 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3405 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3406 "MatConvert_SeqAIJ_SeqCSRPERM", 3407 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3408 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3409 "MatConvert_SeqAIJ_SeqCRL", 3410 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3411 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3412 "MatIsTranspose_SeqAIJ", 3413 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3414 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3415 "MatIsHermitianTranspose_SeqAIJ", 3416 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3417 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3418 "MatSeqAIJSetPreallocation_SeqAIJ", 3419 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3420 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3421 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3422 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3423 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3424 "MatReorderForNonzeroDiagonal_SeqAIJ", 3425 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3426 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3427 "MatMatMult_SeqDense_SeqAIJ", 3428 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3429 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3430 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3431 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3432 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3433 "MatMatMultNumeric_SeqDense_SeqAIJ", 3434 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3435 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3436 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3437 PetscFunctionReturn(0); 3438 } 3439 EXTERN_C_END 3440 3441 #undef __FUNCT__ 3442 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3443 /* 3444 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3445 */ 3446 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3447 { 3448 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3449 PetscErrorCode ierr; 3450 PetscInt i,m = A->rmap->n; 3451 3452 PetscFunctionBegin; 3453 c = (Mat_SeqAIJ*)C->data; 3454 3455 C->factortype = A->factortype; 3456 c->row = 0; 3457 c->col = 0; 3458 c->icol = 0; 3459 c->reallocs = 0; 3460 3461 C->assembled = PETSC_TRUE; 3462 3463 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3464 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3465 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3466 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3467 3468 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3469 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3470 for (i=0; i<m; i++) { 3471 c->imax[i] = a->imax[i]; 3472 c->ilen[i] = a->ilen[i]; 3473 } 3474 3475 /* allocate the matrix space */ 3476 if (mallocmatspace){ 3477 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3478 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3479 c->singlemalloc = PETSC_TRUE; 3480 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3481 if (m > 0) { 3482 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3483 if (cpvalues == MAT_COPY_VALUES) { 3484 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3485 } else { 3486 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3487 } 3488 } 3489 } 3490 3491 c->ignorezeroentries = a->ignorezeroentries; 3492 c->roworiented = a->roworiented; 3493 c->nonew = a->nonew; 3494 if (a->diag) { 3495 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3496 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3497 for (i=0; i<m; i++) { 3498 c->diag[i] = a->diag[i]; 3499 } 3500 } else c->diag = 0; 3501 c->solve_work = 0; 3502 c->saved_values = 0; 3503 c->idiag = 0; 3504 c->ssor_work = 0; 3505 c->keepnonzeropattern = a->keepnonzeropattern; 3506 c->free_a = PETSC_TRUE; 3507 c->free_ij = PETSC_TRUE; 3508 c->xtoy = 0; 3509 c->XtoY = 0; 3510 3511 c->nz = a->nz; 3512 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3513 C->preallocated = PETSC_TRUE; 3514 3515 c->compressedrow.use = a->compressedrow.use; 3516 c->compressedrow.nrows = a->compressedrow.nrows; 3517 c->compressedrow.checked = a->compressedrow.checked; 3518 if (a->compressedrow.checked && a->compressedrow.use){ 3519 i = a->compressedrow.nrows; 3520 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3521 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3522 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3523 } else { 3524 c->compressedrow.use = PETSC_FALSE; 3525 c->compressedrow.i = PETSC_NULL; 3526 c->compressedrow.rindex = PETSC_NULL; 3527 } 3528 C->same_nonzero = A->same_nonzero; 3529 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3530 3531 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3532 PetscFunctionReturn(0); 3533 } 3534 3535 #undef __FUNCT__ 3536 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3537 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3538 { 3539 PetscErrorCode ierr; 3540 3541 PetscFunctionBegin; 3542 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3543 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3544 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3545 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3546 PetscFunctionReturn(0); 3547 } 3548 3549 #undef __FUNCT__ 3550 #define __FUNCT__ "MatLoad_SeqAIJ" 3551 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, const MatType type,Mat *A) 3552 { 3553 Mat_SeqAIJ *a; 3554 Mat B; 3555 PetscErrorCode ierr; 3556 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3557 int fd; 3558 PetscMPIInt size; 3559 MPI_Comm comm; 3560 3561 PetscFunctionBegin; 3562 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3563 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3564 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3565 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3566 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3567 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3568 M = header[1]; N = header[2]; nz = header[3]; 3569 3570 if (nz < 0) { 3571 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3572 } 3573 3574 /* read in row lengths */ 3575 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3576 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3577 3578 /* check if sum of rowlengths is same as nz */ 3579 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3580 if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3581 3582 /* create our matrix */ 3583 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3584 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3585 ierr = MatSetType(B,type);CHKERRQ(ierr); 3586 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3587 a = (Mat_SeqAIJ*)B->data; 3588 3589 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3590 3591 /* read in nonzero values */ 3592 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3593 3594 /* set matrix "i" values */ 3595 a->i[0] = 0; 3596 for (i=1; i<= M; i++) { 3597 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3598 a->ilen[i-1] = rowlengths[i-1]; 3599 } 3600 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3601 3602 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3603 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3604 *A = B; 3605 PetscFunctionReturn(0); 3606 } 3607 3608 #undef __FUNCT__ 3609 #define __FUNCT__ "MatLoadnew_SeqAIJ" 3610 PetscErrorCode MatLoadnew_SeqAIJ(PetscViewer viewer,Mat newMat) 3611 { 3612 Mat_SeqAIJ *a; 3613 PetscErrorCode ierr; 3614 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3615 int fd; 3616 PetscMPIInt size; 3617 MPI_Comm comm; 3618 3619 PetscFunctionBegin; 3620 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3621 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3622 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3623 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3624 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3625 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3626 M = header[1]; N = header[2]; nz = header[3]; 3627 3628 if (nz < 0) { 3629 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3630 } 3631 3632 /* read in row lengths */ 3633 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3634 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3635 3636 /* check if sum of rowlengths is same as nz */ 3637 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3638 if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3639 3640 /* set global size if not set already*/ 3641 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3642 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3643 } else { 3644 /* if sizes and type are already set, check if the vector global sizes are correct */ 3645 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3646 if (M != rows || N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols); 3647 } 3648 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3649 a = (Mat_SeqAIJ*)newMat->data; 3650 3651 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3652 3653 /* read in nonzero values */ 3654 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3655 3656 /* set matrix "i" values */ 3657 a->i[0] = 0; 3658 for (i=1; i<= M; i++) { 3659 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3660 a->ilen[i-1] = rowlengths[i-1]; 3661 } 3662 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3663 3664 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3665 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3666 PetscFunctionReturn(0); 3667 } 3668 3669 #undef __FUNCT__ 3670 #define __FUNCT__ "MatEqual_SeqAIJ" 3671 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3672 { 3673 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3674 PetscErrorCode ierr; 3675 #if defined(PETSC_USE_COMPLEX) 3676 PetscInt k; 3677 #endif 3678 3679 PetscFunctionBegin; 3680 /* If the matrix dimensions are not equal,or no of nonzeros */ 3681 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3682 *flg = PETSC_FALSE; 3683 PetscFunctionReturn(0); 3684 } 3685 3686 /* if the a->i are the same */ 3687 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3688 if (!*flg) PetscFunctionReturn(0); 3689 3690 /* if a->j are the same */ 3691 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3692 if (!*flg) PetscFunctionReturn(0); 3693 3694 /* if a->a are the same */ 3695 #if defined(PETSC_USE_COMPLEX) 3696 for (k=0; k<a->nz; k++){ 3697 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3698 *flg = PETSC_FALSE; 3699 PetscFunctionReturn(0); 3700 } 3701 } 3702 #else 3703 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3704 #endif 3705 PetscFunctionReturn(0); 3706 } 3707 3708 #undef __FUNCT__ 3709 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3710 /*@ 3711 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3712 provided by the user. 3713 3714 Collective on MPI_Comm 3715 3716 Input Parameters: 3717 + comm - must be an MPI communicator of size 1 3718 . m - number of rows 3719 . n - number of columns 3720 . i - row indices 3721 . j - column indices 3722 - a - matrix values 3723 3724 Output Parameter: 3725 . mat - the matrix 3726 3727 Level: intermediate 3728 3729 Notes: 3730 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3731 once the matrix is destroyed 3732 3733 You cannot set new nonzero locations into this matrix, that will generate an error. 3734 3735 The i and j indices are 0 based 3736 3737 The format which is used for the sparse matrix input, is equivalent to a 3738 row-major ordering.. i.e for the following matrix, the input data expected is 3739 as shown: 3740 3741 1 0 0 3742 2 0 3 3743 4 5 6 3744 3745 i = {0,1,3,6} [size = nrow+1 = 3+1] 3746 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3747 v = {1,2,3,4,5,6} [size = nz = 6] 3748 3749 3750 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3751 3752 @*/ 3753 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3754 { 3755 PetscErrorCode ierr; 3756 PetscInt ii; 3757 Mat_SeqAIJ *aij; 3758 #if defined(PETSC_USE_DEBUG) 3759 PetscInt jj; 3760 #endif 3761 3762 PetscFunctionBegin; 3763 if (i[0]) { 3764 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3765 } 3766 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3767 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3768 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3769 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3770 aij = (Mat_SeqAIJ*)(*mat)->data; 3771 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3772 3773 aij->i = i; 3774 aij->j = j; 3775 aij->a = a; 3776 aij->singlemalloc = PETSC_FALSE; 3777 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3778 aij->free_a = PETSC_FALSE; 3779 aij->free_ij = PETSC_FALSE; 3780 3781 for (ii=0; ii<m; ii++) { 3782 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3783 #if defined(PETSC_USE_DEBUG) 3784 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 3785 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3786 if (j[jj] < j[jj-1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is not sorted",jj-i[ii],j[jj],ii); 3787 if (j[jj] == j[jj]-1) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii); 3788 } 3789 #endif 3790 } 3791 #if defined(PETSC_USE_DEBUG) 3792 for (ii=0; ii<aij->i[m]; ii++) { 3793 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3794 if (j[ii] > n - 1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3795 } 3796 #endif 3797 3798 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3799 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3800 PetscFunctionReturn(0); 3801 } 3802 3803 #undef __FUNCT__ 3804 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3805 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3806 { 3807 PetscErrorCode ierr; 3808 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3809 3810 PetscFunctionBegin; 3811 if (coloring->ctype == IS_COLORING_GLOBAL) { 3812 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3813 a->coloring = coloring; 3814 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3815 PetscInt i,*larray; 3816 ISColoring ocoloring; 3817 ISColoringValue *colors; 3818 3819 /* set coloring for diagonal portion */ 3820 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3821 for (i=0; i<A->cmap->n; i++) { 3822 larray[i] = i; 3823 } 3824 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3825 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3826 for (i=0; i<A->cmap->n; i++) { 3827 colors[i] = coloring->colors[larray[i]]; 3828 } 3829 ierr = PetscFree(larray);CHKERRQ(ierr); 3830 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3831 a->coloring = ocoloring; 3832 } 3833 PetscFunctionReturn(0); 3834 } 3835 3836 #if defined(PETSC_HAVE_ADIC) 3837 EXTERN_C_BEGIN 3838 #include "adic/ad_utils.h" 3839 EXTERN_C_END 3840 3841 #undef __FUNCT__ 3842 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3843 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3844 { 3845 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3846 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3847 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3848 ISColoringValue *color; 3849 3850 PetscFunctionBegin; 3851 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3852 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3853 color = a->coloring->colors; 3854 /* loop over rows */ 3855 for (i=0; i<m; i++) { 3856 nz = ii[i+1] - ii[i]; 3857 /* loop over columns putting computed value into matrix */ 3858 for (j=0; j<nz; j++) { 3859 *v++ = values[color[*jj++]]; 3860 } 3861 values += nlen; /* jump to next row of derivatives */ 3862 } 3863 PetscFunctionReturn(0); 3864 } 3865 #endif 3866 3867 #undef __FUNCT__ 3868 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3869 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3870 { 3871 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3872 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3873 MatScalar *v = a->a; 3874 PetscScalar *values = (PetscScalar *)advalues; 3875 ISColoringValue *color; 3876 3877 PetscFunctionBegin; 3878 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3879 color = a->coloring->colors; 3880 /* loop over rows */ 3881 for (i=0; i<m; i++) { 3882 nz = ii[i+1] - ii[i]; 3883 /* loop over columns putting computed value into matrix */ 3884 for (j=0; j<nz; j++) { 3885 *v++ = values[color[*jj++]]; 3886 } 3887 values += nl; /* jump to next row of derivatives */ 3888 } 3889 PetscFunctionReturn(0); 3890 } 3891 3892 /* 3893 Special version for direct calls from Fortran 3894 */ 3895 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3896 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3897 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3898 #define matsetvaluesseqaij_ matsetvaluesseqaij 3899 #endif 3900 3901 /* Change these macros so can be used in void function */ 3902 #undef CHKERRQ 3903 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3904 #undef SETERRQ2 3905 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3906 3907 EXTERN_C_BEGIN 3908 #undef __FUNCT__ 3909 #define __FUNCT__ "matsetvaluesseqaij_" 3910 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3911 { 3912 Mat A = *AA; 3913 PetscInt m = *mm, n = *nn; 3914 InsertMode is = *isis; 3915 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3916 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3917 PetscInt *imax,*ai,*ailen; 3918 PetscErrorCode ierr; 3919 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3920 MatScalar *ap,value,*aa; 3921 PetscTruth ignorezeroentries = a->ignorezeroentries; 3922 PetscTruth roworiented = a->roworiented; 3923 3924 PetscFunctionBegin; 3925 ierr = MatPreallocated(A);CHKERRQ(ierr); 3926 imax = a->imax; 3927 ai = a->i; 3928 ailen = a->ilen; 3929 aj = a->j; 3930 aa = a->a; 3931 3932 for (k=0; k<m; k++) { /* loop over added rows */ 3933 row = im[k]; 3934 if (row < 0) continue; 3935 #if defined(PETSC_USE_DEBUG) 3936 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3937 #endif 3938 rp = aj + ai[row]; ap = aa + ai[row]; 3939 rmax = imax[row]; nrow = ailen[row]; 3940 low = 0; 3941 high = nrow; 3942 for (l=0; l<n; l++) { /* loop over added columns */ 3943 if (in[l] < 0) continue; 3944 #if defined(PETSC_USE_DEBUG) 3945 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3946 #endif 3947 col = in[l]; 3948 if (roworiented) { 3949 value = v[l + k*n]; 3950 } else { 3951 value = v[k + l*m]; 3952 } 3953 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3954 3955 if (col <= lastcol) low = 0; else high = nrow; 3956 lastcol = col; 3957 while (high-low > 5) { 3958 t = (low+high)/2; 3959 if (rp[t] > col) high = t; 3960 else low = t; 3961 } 3962 for (i=low; i<high; i++) { 3963 if (rp[i] > col) break; 3964 if (rp[i] == col) { 3965 if (is == ADD_VALUES) ap[i] += value; 3966 else ap[i] = value; 3967 goto noinsert; 3968 } 3969 } 3970 if (value == 0.0 && ignorezeroentries) goto noinsert; 3971 if (nonew == 1) goto noinsert; 3972 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3973 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3974 N = nrow++ - 1; a->nz++; high++; 3975 /* shift up all the later entries in this row */ 3976 for (ii=N; ii>=i; ii--) { 3977 rp[ii+1] = rp[ii]; 3978 ap[ii+1] = ap[ii]; 3979 } 3980 rp[i] = col; 3981 ap[i] = value; 3982 noinsert:; 3983 low = i + 1; 3984 } 3985 ailen[row] = nrow; 3986 } 3987 A->same_nonzero = PETSC_FALSE; 3988 PetscFunctionReturnVoid(); 3989 } 3990 EXTERN_C_END 3991