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 = VecCUDAAllocateCheck(yy);CHKERRQ(ierr); 1031 cusp::multiply(a->GPUmatrix,xx->GPUarray,yy->GPUarray); 1032 yy->valid_GPU_array = PETSC_CUDA_GPU; 1033 #else 1034 for (i=0; i<m; i++) { 1035 n = ii[i+1] - ii[i]; 1036 aj = a->j + ii[i]; 1037 aa = a->a + ii[i]; 1038 sum = 0.0; 1039 nonzerorow += (n>0); 1040 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1041 y[i] = sum; 1042 } 1043 #endif 1044 #endif 1045 } 1046 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1047 #if !defined(PETSC_HAVE_CUDA) 1048 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1049 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1050 #endif 1051 PetscFunctionReturn(0); 1052 } 1053 1054 #include "../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h" 1055 #undef __FUNCT__ 1056 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1057 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1058 { 1059 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1060 PetscScalar *x,*y,*z; 1061 const MatScalar *aa; 1062 PetscErrorCode ierr; 1063 PetscInt m = A->rmap->n,*aj,*ii; 1064 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1065 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 1066 PetscScalar sum; 1067 PetscTruth usecprow=a->compressedrow.use; 1068 #endif 1069 1070 PetscFunctionBegin; 1071 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1072 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1073 if (zz != yy) { 1074 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1075 } else { 1076 z = y; 1077 } 1078 1079 aj = a->j; 1080 aa = a->a; 1081 ii = a->i; 1082 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1083 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1084 #else 1085 if (usecprow){ /* use compressed row format */ 1086 if (zz != yy){ 1087 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1088 } 1089 m = a->compressedrow.nrows; 1090 ii = a->compressedrow.i; 1091 ridx = a->compressedrow.rindex; 1092 for (i=0; i<m; i++){ 1093 n = ii[i+1] - ii[i]; 1094 aj = a->j + ii[i]; 1095 aa = a->a + ii[i]; 1096 sum = y[*ridx]; 1097 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 1098 z[*ridx++] = sum; 1099 } 1100 } else { /* do not use compressed row format */ 1101 for (i=0; i<m; i++) { 1102 jrow = ii[i]; 1103 n = ii[i+1] - jrow; 1104 sum = y[i]; 1105 for (j=0; j<n; j++) { 1106 sum += aa[jrow]*x[aj[jrow]]; jrow++; 1107 } 1108 z[i] = sum; 1109 } 1110 } 1111 #endif 1112 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1113 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1114 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1115 if (zz != yy) { 1116 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1117 } 1118 PetscFunctionReturn(0); 1119 } 1120 1121 /* 1122 Adds diagonal pointers to sparse matrix structure. 1123 */ 1124 #undef __FUNCT__ 1125 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1126 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1127 { 1128 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1129 PetscErrorCode ierr; 1130 PetscInt i,j,m = A->rmap->n; 1131 1132 PetscFunctionBegin; 1133 if (!a->diag) { 1134 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1135 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1136 } 1137 for (i=0; i<A->rmap->n; i++) { 1138 a->diag[i] = a->i[i+1]; 1139 for (j=a->i[i]; j<a->i[i+1]; j++) { 1140 if (a->j[j] == i) { 1141 a->diag[i] = j; 1142 break; 1143 } 1144 } 1145 } 1146 PetscFunctionReturn(0); 1147 } 1148 1149 /* 1150 Checks for missing diagonals 1151 */ 1152 #undef __FUNCT__ 1153 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1154 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1155 { 1156 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1157 PetscInt *diag,*jj = a->j,i; 1158 1159 PetscFunctionBegin; 1160 *missing = PETSC_FALSE; 1161 if (A->rmap->n > 0 && !jj) { 1162 *missing = PETSC_TRUE; 1163 if (d) *d = 0; 1164 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1165 } else { 1166 diag = a->diag; 1167 for (i=0; i<A->rmap->n; i++) { 1168 if (jj[diag[i]] != i) { 1169 *missing = PETSC_TRUE; 1170 if (d) *d = i; 1171 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1172 } 1173 } 1174 } 1175 PetscFunctionReturn(0); 1176 } 1177 1178 EXTERN_C_BEGIN 1179 #undef __FUNCT__ 1180 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1181 PetscErrorCode PETSCMAT_DLLEXPORT MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1182 { 1183 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1184 PetscErrorCode ierr; 1185 PetscInt i,*diag,m = A->rmap->n; 1186 MatScalar *v = a->a; 1187 PetscScalar *idiag,*mdiag; 1188 1189 PetscFunctionBegin; 1190 if (a->idiagvalid) PetscFunctionReturn(0); 1191 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1192 diag = a->diag; 1193 if (!a->idiag) { 1194 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1195 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1196 v = a->a; 1197 } 1198 mdiag = a->mdiag; 1199 idiag = a->idiag; 1200 1201 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1202 for (i=0; i<m; i++) { 1203 mdiag[i] = v[diag[i]]; 1204 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1205 idiag[i] = 1.0/v[diag[i]]; 1206 } 1207 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1208 } else { 1209 for (i=0; i<m; i++) { 1210 mdiag[i] = v[diag[i]]; 1211 idiag[i] = omega/(fshift + v[diag[i]]); 1212 } 1213 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1214 } 1215 a->idiagvalid = PETSC_TRUE; 1216 PetscFunctionReturn(0); 1217 } 1218 EXTERN_C_END 1219 1220 #include "../src/mat/impls/aij/seq/ftn-kernels/frelax.h" 1221 #undef __FUNCT__ 1222 #define __FUNCT__ "MatSOR_SeqAIJ" 1223 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1224 { 1225 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1226 PetscScalar *x,d,sum,*t,scale; 1227 const MatScalar *v = a->a,*idiag=0,*mdiag; 1228 const PetscScalar *b, *bs,*xb, *ts; 1229 PetscErrorCode ierr; 1230 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1231 const PetscInt *idx,*diag; 1232 1233 PetscFunctionBegin; 1234 its = its*lits; 1235 1236 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1237 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1238 a->fshift = fshift; 1239 a->omega = omega; 1240 1241 diag = a->diag; 1242 t = a->ssor_work; 1243 idiag = a->idiag; 1244 mdiag = a->mdiag; 1245 1246 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1247 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1248 CHKMEMQ; 1249 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1250 if (flag == SOR_APPLY_UPPER) { 1251 /* apply (U + D/omega) to the vector */ 1252 bs = b; 1253 for (i=0; i<m; i++) { 1254 d = fshift + mdiag[i]; 1255 n = a->i[i+1] - diag[i] - 1; 1256 idx = a->j + diag[i] + 1; 1257 v = a->a + diag[i] + 1; 1258 sum = b[i]*d/omega; 1259 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1260 x[i] = sum; 1261 } 1262 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1263 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1264 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1265 PetscFunctionReturn(0); 1266 } 1267 1268 if (flag == SOR_APPLY_LOWER) { 1269 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1270 } else if (flag & SOR_EISENSTAT) { 1271 /* Let A = L + U + D; where L is lower trianglar, 1272 U is upper triangular, E = D/omega; This routine applies 1273 1274 (L + E)^{-1} A (U + E)^{-1} 1275 1276 to a vector efficiently using Eisenstat's trick. 1277 */ 1278 scale = (2.0/omega) - 1.0; 1279 1280 /* x = (E + U)^{-1} b */ 1281 for (i=m-1; i>=0; i--) { 1282 n = a->i[i+1] - diag[i] - 1; 1283 idx = a->j + diag[i] + 1; 1284 v = a->a + diag[i] + 1; 1285 sum = b[i]; 1286 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1287 x[i] = sum*idiag[i]; 1288 } 1289 1290 /* t = b - (2*E - D)x */ 1291 v = a->a; 1292 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1293 1294 /* t = (E + L)^{-1}t */ 1295 ts = t; 1296 diag = a->diag; 1297 for (i=0; i<m; i++) { 1298 n = diag[i] - a->i[i]; 1299 idx = a->j + a->i[i]; 1300 v = a->a + a->i[i]; 1301 sum = t[i]; 1302 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1303 t[i] = sum*idiag[i]; 1304 /* x = x + t */ 1305 x[i] += t[i]; 1306 } 1307 1308 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1309 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1310 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1311 PetscFunctionReturn(0); 1312 } 1313 if (flag & SOR_ZERO_INITIAL_GUESS) { 1314 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1315 for (i=0; i<m; i++) { 1316 n = diag[i] - a->i[i]; 1317 idx = a->j + a->i[i]; 1318 v = a->a + a->i[i]; 1319 sum = b[i]; 1320 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1321 t[i] = sum; 1322 x[i] = sum*idiag[i]; 1323 } 1324 xb = t; 1325 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1326 } else xb = b; 1327 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1328 for (i=m-1; i>=0; i--) { 1329 n = a->i[i+1] - diag[i] - 1; 1330 idx = a->j + diag[i] + 1; 1331 v = a->a + diag[i] + 1; 1332 sum = xb[i]; 1333 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1334 if (xb == b) { 1335 x[i] = sum*idiag[i]; 1336 } else { 1337 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1338 } 1339 } 1340 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1341 } 1342 its--; 1343 } 1344 while (its--) { 1345 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1346 for (i=0; i<m; i++) { 1347 n = a->i[i+1] - a->i[i]; 1348 idx = a->j + a->i[i]; 1349 v = a->a + a->i[i]; 1350 sum = b[i]; 1351 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1352 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1353 } 1354 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1355 } 1356 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1357 for (i=m-1; i>=0; i--) { 1358 n = a->i[i+1] - a->i[i]; 1359 idx = a->j + a->i[i]; 1360 v = a->a + a->i[i]; 1361 sum = b[i]; 1362 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1363 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1364 } 1365 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1366 } 1367 } 1368 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1369 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1370 CHKMEMQ; PetscFunctionReturn(0); 1371 } 1372 1373 1374 #undef __FUNCT__ 1375 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1376 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1377 { 1378 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1379 1380 PetscFunctionBegin; 1381 info->block_size = 1.0; 1382 info->nz_allocated = (double)a->maxnz; 1383 info->nz_used = (double)a->nz; 1384 info->nz_unneeded = (double)(a->maxnz - a->nz); 1385 info->assemblies = (double)A->num_ass; 1386 info->mallocs = (double)A->info.mallocs; 1387 info->memory = ((PetscObject)A)->mem; 1388 if (A->factortype) { 1389 info->fill_ratio_given = A->info.fill_ratio_given; 1390 info->fill_ratio_needed = A->info.fill_ratio_needed; 1391 info->factor_mallocs = A->info.factor_mallocs; 1392 } else { 1393 info->fill_ratio_given = 0; 1394 info->fill_ratio_needed = 0; 1395 info->factor_mallocs = 0; 1396 } 1397 PetscFunctionReturn(0); 1398 } 1399 1400 #undef __FUNCT__ 1401 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1402 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1403 { 1404 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1405 PetscInt i,m = A->rmap->n - 1,d = 0; 1406 PetscErrorCode ierr; 1407 PetscTruth missing; 1408 1409 PetscFunctionBegin; 1410 if (a->keepnonzeropattern) { 1411 for (i=0; i<N; i++) { 1412 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1413 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1414 } 1415 if (diag != 0.0) { 1416 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1417 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1418 for (i=0; i<N; i++) { 1419 a->a[a->diag[rows[i]]] = diag; 1420 } 1421 } 1422 A->same_nonzero = PETSC_TRUE; 1423 } else { 1424 if (diag != 0.0) { 1425 for (i=0; i<N; i++) { 1426 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1427 if (a->ilen[rows[i]] > 0) { 1428 a->ilen[rows[i]] = 1; 1429 a->a[a->i[rows[i]]] = diag; 1430 a->j[a->i[rows[i]]] = rows[i]; 1431 } else { /* in case row was completely empty */ 1432 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1433 } 1434 } 1435 } else { 1436 for (i=0; i<N; i++) { 1437 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1438 a->ilen[rows[i]] = 0; 1439 } 1440 } 1441 A->same_nonzero = PETSC_FALSE; 1442 } 1443 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1444 PetscFunctionReturn(0); 1445 } 1446 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "MatGetRow_SeqAIJ" 1449 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1450 { 1451 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1452 PetscInt *itmp; 1453 1454 PetscFunctionBegin; 1455 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1456 1457 *nz = a->i[row+1] - a->i[row]; 1458 if (v) *v = a->a + a->i[row]; 1459 if (idx) { 1460 itmp = a->j + a->i[row]; 1461 if (*nz) { 1462 *idx = itmp; 1463 } 1464 else *idx = 0; 1465 } 1466 PetscFunctionReturn(0); 1467 } 1468 1469 /* remove this function? */ 1470 #undef __FUNCT__ 1471 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1472 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1473 { 1474 PetscFunctionBegin; 1475 PetscFunctionReturn(0); 1476 } 1477 1478 #undef __FUNCT__ 1479 #define __FUNCT__ "MatNorm_SeqAIJ" 1480 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1481 { 1482 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1483 MatScalar *v = a->a; 1484 PetscReal sum = 0.0; 1485 PetscErrorCode ierr; 1486 PetscInt i,j; 1487 1488 PetscFunctionBegin; 1489 if (type == NORM_FROBENIUS) { 1490 for (i=0; i<a->nz; i++) { 1491 #if defined(PETSC_USE_COMPLEX) 1492 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1493 #else 1494 sum += (*v)*(*v); v++; 1495 #endif 1496 } 1497 *nrm = sqrt(sum); 1498 } else if (type == NORM_1) { 1499 PetscReal *tmp; 1500 PetscInt *jj = a->j; 1501 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1502 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1503 *nrm = 0.0; 1504 for (j=0; j<a->nz; j++) { 1505 tmp[*jj++] += PetscAbsScalar(*v); v++; 1506 } 1507 for (j=0; j<A->cmap->n; j++) { 1508 if (tmp[j] > *nrm) *nrm = tmp[j]; 1509 } 1510 ierr = PetscFree(tmp);CHKERRQ(ierr); 1511 } else if (type == NORM_INFINITY) { 1512 *nrm = 0.0; 1513 for (j=0; j<A->rmap->n; j++) { 1514 v = a->a + a->i[j]; 1515 sum = 0.0; 1516 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1517 sum += PetscAbsScalar(*v); v++; 1518 } 1519 if (sum > *nrm) *nrm = sum; 1520 } 1521 } else { 1522 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1523 } 1524 PetscFunctionReturn(0); 1525 } 1526 1527 #undef __FUNCT__ 1528 #define __FUNCT__ "MatTranspose_SeqAIJ" 1529 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1530 { 1531 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1532 Mat C; 1533 PetscErrorCode ierr; 1534 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1535 MatScalar *array = a->a; 1536 1537 PetscFunctionBegin; 1538 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"); 1539 1540 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1541 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1542 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1543 1544 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1545 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1546 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1547 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1548 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1549 ierr = PetscFree(col);CHKERRQ(ierr); 1550 } else { 1551 C = *B; 1552 } 1553 1554 for (i=0; i<m; i++) { 1555 len = ai[i+1]-ai[i]; 1556 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1557 array += len; 1558 aj += len; 1559 } 1560 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1561 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1562 1563 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1564 *B = C; 1565 } else { 1566 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1567 } 1568 PetscFunctionReturn(0); 1569 } 1570 1571 EXTERN_C_BEGIN 1572 #undef __FUNCT__ 1573 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1574 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1575 { 1576 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1577 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1578 MatScalar *va,*vb; 1579 PetscErrorCode ierr; 1580 PetscInt ma,na,mb,nb, i; 1581 1582 PetscFunctionBegin; 1583 bij = (Mat_SeqAIJ *) B->data; 1584 1585 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1586 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1587 if (ma!=nb || na!=mb){ 1588 *f = PETSC_FALSE; 1589 PetscFunctionReturn(0); 1590 } 1591 aii = aij->i; bii = bij->i; 1592 adx = aij->j; bdx = bij->j; 1593 va = aij->a; vb = bij->a; 1594 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1595 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1596 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1597 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1598 1599 *f = PETSC_TRUE; 1600 for (i=0; i<ma; i++) { 1601 while (aptr[i]<aii[i+1]) { 1602 PetscInt idc,idr; 1603 PetscScalar vc,vr; 1604 /* column/row index/value */ 1605 idc = adx[aptr[i]]; 1606 idr = bdx[bptr[idc]]; 1607 vc = va[aptr[i]]; 1608 vr = vb[bptr[idc]]; 1609 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1610 *f = PETSC_FALSE; 1611 goto done; 1612 } else { 1613 aptr[i]++; 1614 if (B || i!=idc) bptr[idc]++; 1615 } 1616 } 1617 } 1618 done: 1619 ierr = PetscFree(aptr);CHKERRQ(ierr); 1620 if (B) { 1621 ierr = PetscFree(bptr);CHKERRQ(ierr); 1622 } 1623 PetscFunctionReturn(0); 1624 } 1625 EXTERN_C_END 1626 1627 EXTERN_C_BEGIN 1628 #undef __FUNCT__ 1629 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1630 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1631 { 1632 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1633 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1634 MatScalar *va,*vb; 1635 PetscErrorCode ierr; 1636 PetscInt ma,na,mb,nb, i; 1637 1638 PetscFunctionBegin; 1639 bij = (Mat_SeqAIJ *) B->data; 1640 1641 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1642 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1643 if (ma!=nb || na!=mb){ 1644 *f = PETSC_FALSE; 1645 PetscFunctionReturn(0); 1646 } 1647 aii = aij->i; bii = bij->i; 1648 adx = aij->j; bdx = bij->j; 1649 va = aij->a; vb = bij->a; 1650 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1651 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1652 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1653 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1654 1655 *f = PETSC_TRUE; 1656 for (i=0; i<ma; i++) { 1657 while (aptr[i]<aii[i+1]) { 1658 PetscInt idc,idr; 1659 PetscScalar vc,vr; 1660 /* column/row index/value */ 1661 idc = adx[aptr[i]]; 1662 idr = bdx[bptr[idc]]; 1663 vc = va[aptr[i]]; 1664 vr = vb[bptr[idc]]; 1665 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1666 *f = PETSC_FALSE; 1667 goto done; 1668 } else { 1669 aptr[i]++; 1670 if (B || i!=idc) bptr[idc]++; 1671 } 1672 } 1673 } 1674 done: 1675 ierr = PetscFree(aptr);CHKERRQ(ierr); 1676 if (B) { 1677 ierr = PetscFree(bptr);CHKERRQ(ierr); 1678 } 1679 PetscFunctionReturn(0); 1680 } 1681 EXTERN_C_END 1682 1683 #undef __FUNCT__ 1684 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1685 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1686 { 1687 PetscErrorCode ierr; 1688 PetscFunctionBegin; 1689 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1690 PetscFunctionReturn(0); 1691 } 1692 1693 #undef __FUNCT__ 1694 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1695 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1696 { 1697 PetscErrorCode ierr; 1698 PetscFunctionBegin; 1699 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1700 PetscFunctionReturn(0); 1701 } 1702 1703 #undef __FUNCT__ 1704 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1705 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1706 { 1707 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1708 PetscScalar *l,*r,x; 1709 MatScalar *v; 1710 PetscErrorCode ierr; 1711 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1712 1713 PetscFunctionBegin; 1714 if (ll) { 1715 /* The local size is used so that VecMPI can be passed to this routine 1716 by MatDiagonalScale_MPIAIJ */ 1717 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1718 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1719 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1720 v = a->a; 1721 for (i=0; i<m; i++) { 1722 x = l[i]; 1723 M = a->i[i+1] - a->i[i]; 1724 for (j=0; j<M; j++) { (*v++) *= x;} 1725 } 1726 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1727 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1728 } 1729 if (rr) { 1730 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1731 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1732 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1733 v = a->a; jj = a->j; 1734 for (i=0; i<nz; i++) { 1735 (*v++) *= r[*jj++]; 1736 } 1737 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1738 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1739 } 1740 PetscFunctionReturn(0); 1741 } 1742 1743 #undef __FUNCT__ 1744 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1745 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1746 { 1747 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1748 PetscErrorCode ierr; 1749 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1750 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1751 const PetscInt *irow,*icol; 1752 PetscInt nrows,ncols; 1753 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1754 MatScalar *a_new,*mat_a; 1755 Mat C; 1756 PetscTruth stride,sorted; 1757 1758 PetscFunctionBegin; 1759 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1760 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1761 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1762 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1763 1764 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1765 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1766 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1767 1768 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1769 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1770 if (stride && step == 1) { 1771 /* special case of contiguous rows */ 1772 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1773 /* loop over new rows determining lens and starting points */ 1774 for (i=0; i<nrows; i++) { 1775 kstart = ai[irow[i]]; 1776 kend = kstart + ailen[irow[i]]; 1777 for (k=kstart; k<kend; k++) { 1778 if (aj[k] >= first) { 1779 starts[i] = k; 1780 break; 1781 } 1782 } 1783 sum = 0; 1784 while (k < kend) { 1785 if (aj[k++] >= first+ncols) break; 1786 sum++; 1787 } 1788 lens[i] = sum; 1789 } 1790 /* create submatrix */ 1791 if (scall == MAT_REUSE_MATRIX) { 1792 PetscInt n_cols,n_rows; 1793 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1794 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1795 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1796 C = *B; 1797 } else { 1798 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1799 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1800 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1801 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1802 } 1803 c = (Mat_SeqAIJ*)C->data; 1804 1805 /* loop over rows inserting into submatrix */ 1806 a_new = c->a; 1807 j_new = c->j; 1808 i_new = c->i; 1809 1810 for (i=0; i<nrows; i++) { 1811 ii = starts[i]; 1812 lensi = lens[i]; 1813 for (k=0; k<lensi; k++) { 1814 *j_new++ = aj[ii+k] - first; 1815 } 1816 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1817 a_new += lensi; 1818 i_new[i+1] = i_new[i] + lensi; 1819 c->ilen[i] = lensi; 1820 } 1821 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1822 } else { 1823 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1824 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1825 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1826 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1827 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1828 /* determine lens of each row */ 1829 for (i=0; i<nrows; i++) { 1830 kstart = ai[irow[i]]; 1831 kend = kstart + a->ilen[irow[i]]; 1832 lens[i] = 0; 1833 for (k=kstart; k<kend; k++) { 1834 if (smap[aj[k]]) { 1835 lens[i]++; 1836 } 1837 } 1838 } 1839 /* Create and fill new matrix */ 1840 if (scall == MAT_REUSE_MATRIX) { 1841 PetscTruth equal; 1842 1843 c = (Mat_SeqAIJ *)((*B)->data); 1844 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1845 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1846 if (!equal) { 1847 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1848 } 1849 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1850 C = *B; 1851 } else { 1852 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1853 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1854 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1855 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1856 } 1857 c = (Mat_SeqAIJ *)(C->data); 1858 for (i=0; i<nrows; i++) { 1859 row = irow[i]; 1860 kstart = ai[row]; 1861 kend = kstart + a->ilen[row]; 1862 mat_i = c->i[i]; 1863 mat_j = c->j + mat_i; 1864 mat_a = c->a + mat_i; 1865 mat_ilen = c->ilen + i; 1866 for (k=kstart; k<kend; k++) { 1867 if ((tcol=smap[a->j[k]])) { 1868 *mat_j++ = tcol - 1; 1869 *mat_a++ = a->a[k]; 1870 (*mat_ilen)++; 1871 1872 } 1873 } 1874 } 1875 /* Free work space */ 1876 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1877 ierr = PetscFree(smap);CHKERRQ(ierr); 1878 ierr = PetscFree(lens);CHKERRQ(ierr); 1879 } 1880 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1881 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1882 1883 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1884 *B = C; 1885 PetscFunctionReturn(0); 1886 } 1887 1888 #undef __FUNCT__ 1889 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1890 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1891 { 1892 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1893 PetscErrorCode ierr; 1894 Mat outA; 1895 PetscTruth row_identity,col_identity; 1896 1897 PetscFunctionBegin; 1898 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1899 1900 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1901 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1902 1903 outA = inA; 1904 outA->factortype = MAT_FACTOR_LU; 1905 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1906 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1907 a->row = row; 1908 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1909 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1910 a->col = col; 1911 1912 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1913 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1914 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1915 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1916 1917 if (!a->solve_work) { /* this matrix may have been factored before */ 1918 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1919 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1920 } 1921 1922 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1923 if (row_identity && col_identity) { 1924 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 1925 } else { 1926 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1927 } 1928 PetscFunctionReturn(0); 1929 } 1930 1931 #undef __FUNCT__ 1932 #define __FUNCT__ "MatScale_SeqAIJ" 1933 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1934 { 1935 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1936 PetscScalar oalpha = alpha; 1937 PetscErrorCode ierr; 1938 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1939 1940 PetscFunctionBegin; 1941 BLASscal_(&bnz,&oalpha,a->a,&one); 1942 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1943 PetscFunctionReturn(0); 1944 } 1945 1946 #undef __FUNCT__ 1947 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1948 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1949 { 1950 PetscErrorCode ierr; 1951 PetscInt i; 1952 1953 PetscFunctionBegin; 1954 if (scall == MAT_INITIAL_MATRIX) { 1955 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1956 } 1957 1958 for (i=0; i<n; i++) { 1959 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1960 } 1961 PetscFunctionReturn(0); 1962 } 1963 1964 #undef __FUNCT__ 1965 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1966 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1967 { 1968 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1969 PetscErrorCode ierr; 1970 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1971 const PetscInt *idx; 1972 PetscInt start,end,*ai,*aj; 1973 PetscBT table; 1974 1975 PetscFunctionBegin; 1976 m = A->rmap->n; 1977 ai = a->i; 1978 aj = a->j; 1979 1980 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1981 1982 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1983 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1984 1985 for (i=0; i<is_max; i++) { 1986 /* Initialize the two local arrays */ 1987 isz = 0; 1988 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1989 1990 /* Extract the indices, assume there can be duplicate entries */ 1991 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1992 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1993 1994 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1995 for (j=0; j<n ; ++j){ 1996 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1997 } 1998 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1999 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 2000 2001 k = 0; 2002 for (j=0; j<ov; j++){ /* for each overlap */ 2003 n = isz; 2004 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 2005 row = nidx[k]; 2006 start = ai[row]; 2007 end = ai[row+1]; 2008 for (l = start; l<end ; l++){ 2009 val = aj[l] ; 2010 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 2011 } 2012 } 2013 } 2014 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 2015 } 2016 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 2017 ierr = PetscFree(nidx);CHKERRQ(ierr); 2018 PetscFunctionReturn(0); 2019 } 2020 2021 /* -------------------------------------------------------------- */ 2022 #undef __FUNCT__ 2023 #define __FUNCT__ "MatPermute_SeqAIJ" 2024 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2025 { 2026 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2027 PetscErrorCode ierr; 2028 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2029 const PetscInt *row,*col; 2030 PetscInt *cnew,j,*lens; 2031 IS icolp,irowp; 2032 PetscInt *cwork = PETSC_NULL; 2033 PetscScalar *vwork = PETSC_NULL; 2034 2035 PetscFunctionBegin; 2036 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2037 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2038 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2039 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2040 2041 /* determine lengths of permuted rows */ 2042 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2043 for (i=0; i<m; i++) { 2044 lens[row[i]] = a->i[i+1] - a->i[i]; 2045 } 2046 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2047 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2048 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2049 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2050 ierr = PetscFree(lens);CHKERRQ(ierr); 2051 2052 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2053 for (i=0; i<m; i++) { 2054 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2055 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2056 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2057 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2058 } 2059 ierr = PetscFree(cnew);CHKERRQ(ierr); 2060 (*B)->assembled = PETSC_FALSE; 2061 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2062 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2063 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2064 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2065 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2066 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2067 PetscFunctionReturn(0); 2068 } 2069 2070 #undef __FUNCT__ 2071 #define __FUNCT__ "MatCopy_SeqAIJ" 2072 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2073 { 2074 PetscErrorCode ierr; 2075 2076 PetscFunctionBegin; 2077 /* If the two matrices have the same copy implementation, use fast copy. */ 2078 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2079 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2080 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2081 2082 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2083 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2084 } 2085 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2086 } else { 2087 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2088 } 2089 PetscFunctionReturn(0); 2090 } 2091 2092 #undef __FUNCT__ 2093 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2094 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2095 { 2096 PetscErrorCode ierr; 2097 2098 PetscFunctionBegin; 2099 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2100 PetscFunctionReturn(0); 2101 } 2102 2103 #undef __FUNCT__ 2104 #define __FUNCT__ "MatGetArray_SeqAIJ" 2105 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2106 { 2107 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2108 PetscFunctionBegin; 2109 *array = a->a; 2110 PetscFunctionReturn(0); 2111 } 2112 2113 #undef __FUNCT__ 2114 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2115 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2116 { 2117 PetscFunctionBegin; 2118 PetscFunctionReturn(0); 2119 } 2120 2121 #undef __FUNCT__ 2122 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2123 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2124 { 2125 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2126 PetscErrorCode ierr; 2127 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2128 PetscScalar dx,*y,*xx,*w3_array; 2129 PetscScalar *vscale_array; 2130 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2131 Vec w1,w2,w3; 2132 void *fctx = coloring->fctx; 2133 PetscTruth flg = PETSC_FALSE; 2134 2135 PetscFunctionBegin; 2136 if (!coloring->w1) { 2137 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2138 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2139 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2140 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2141 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2142 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2143 } 2144 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2145 2146 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2147 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2148 if (flg) { 2149 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2150 } else { 2151 PetscTruth assembled; 2152 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2153 if (assembled) { 2154 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2155 } 2156 } 2157 2158 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2159 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2160 2161 /* 2162 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2163 coloring->F for the coarser grids from the finest 2164 */ 2165 if (coloring->F) { 2166 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2167 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2168 if (m1 != m2) { 2169 coloring->F = 0; 2170 } 2171 } 2172 2173 if (coloring->F) { 2174 w1 = coloring->F; 2175 coloring->F = 0; 2176 } else { 2177 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2178 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2179 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2180 } 2181 2182 /* 2183 Compute all the scale factors and share with other processors 2184 */ 2185 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2186 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2187 for (k=0; k<coloring->ncolors; k++) { 2188 /* 2189 Loop over each column associated with color adding the 2190 perturbation to the vector w3. 2191 */ 2192 for (l=0; l<coloring->ncolumns[k]; l++) { 2193 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2194 dx = xx[col]; 2195 if (dx == 0.0) dx = 1.0; 2196 #if !defined(PETSC_USE_COMPLEX) 2197 if (dx < umin && dx >= 0.0) dx = umin; 2198 else if (dx < 0.0 && dx > -umin) dx = -umin; 2199 #else 2200 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2201 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2202 #endif 2203 dx *= epsilon; 2204 vscale_array[col] = 1.0/dx; 2205 } 2206 } 2207 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2208 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2209 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2210 2211 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2212 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2213 2214 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2215 else vscaleforrow = coloring->columnsforrow; 2216 2217 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2218 /* 2219 Loop over each color 2220 */ 2221 for (k=0; k<coloring->ncolors; k++) { 2222 coloring->currentcolor = k; 2223 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2224 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2225 /* 2226 Loop over each column associated with color adding the 2227 perturbation to the vector w3. 2228 */ 2229 for (l=0; l<coloring->ncolumns[k]; l++) { 2230 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2231 dx = xx[col]; 2232 if (dx == 0.0) dx = 1.0; 2233 #if !defined(PETSC_USE_COMPLEX) 2234 if (dx < umin && dx >= 0.0) dx = umin; 2235 else if (dx < 0.0 && dx > -umin) dx = -umin; 2236 #else 2237 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2238 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2239 #endif 2240 dx *= epsilon; 2241 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2242 w3_array[col] += dx; 2243 } 2244 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2245 2246 /* 2247 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2248 */ 2249 2250 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2251 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2252 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2253 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2254 2255 /* 2256 Loop over rows of vector, putting results into Jacobian matrix 2257 */ 2258 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2259 for (l=0; l<coloring->nrows[k]; l++) { 2260 row = coloring->rows[k][l]; 2261 col = coloring->columnsforrow[k][l]; 2262 y[row] *= vscale_array[vscaleforrow[k][l]]; 2263 srow = row + start; 2264 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2265 } 2266 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2267 } 2268 coloring->currentcolor = k; 2269 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2270 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2271 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2272 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2273 PetscFunctionReturn(0); 2274 } 2275 2276 #undef __FUNCT__ 2277 #define __FUNCT__ "MatAXPY_SeqAIJ" 2278 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2279 { 2280 PetscErrorCode ierr; 2281 PetscInt i; 2282 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2283 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2284 2285 PetscFunctionBegin; 2286 if (str == SAME_NONZERO_PATTERN) { 2287 PetscScalar alpha = a; 2288 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2289 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2290 if (y->xtoy && y->XtoY != X) { 2291 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2292 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2293 } 2294 if (!y->xtoy) { /* get xtoy */ 2295 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2296 y->XtoY = X; 2297 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2298 } 2299 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2300 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2301 } else { 2302 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2303 } 2304 PetscFunctionReturn(0); 2305 } 2306 2307 #undef __FUNCT__ 2308 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2309 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2310 { 2311 PetscErrorCode ierr; 2312 2313 PetscFunctionBegin; 2314 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2315 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2316 PetscFunctionReturn(0); 2317 } 2318 2319 #undef __FUNCT__ 2320 #define __FUNCT__ "MatConjugate_SeqAIJ" 2321 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2322 { 2323 #if defined(PETSC_USE_COMPLEX) 2324 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2325 PetscInt i,nz; 2326 PetscScalar *a; 2327 2328 PetscFunctionBegin; 2329 nz = aij->nz; 2330 a = aij->a; 2331 for (i=0; i<nz; i++) { 2332 a[i] = PetscConj(a[i]); 2333 } 2334 #else 2335 PetscFunctionBegin; 2336 #endif 2337 PetscFunctionReturn(0); 2338 } 2339 2340 #undef __FUNCT__ 2341 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2342 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2343 { 2344 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2345 PetscErrorCode ierr; 2346 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2347 PetscReal atmp; 2348 PetscScalar *x; 2349 MatScalar *aa; 2350 2351 PetscFunctionBegin; 2352 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2353 aa = a->a; 2354 ai = a->i; 2355 aj = a->j; 2356 2357 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2358 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2359 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2360 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2361 for (i=0; i<m; i++) { 2362 ncols = ai[1] - ai[0]; ai++; 2363 x[i] = 0.0; 2364 for (j=0; j<ncols; j++){ 2365 atmp = PetscAbsScalar(*aa); 2366 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2367 aa++; aj++; 2368 } 2369 } 2370 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2371 PetscFunctionReturn(0); 2372 } 2373 2374 #undef __FUNCT__ 2375 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2376 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2377 { 2378 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2379 PetscErrorCode ierr; 2380 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2381 PetscScalar *x; 2382 MatScalar *aa; 2383 2384 PetscFunctionBegin; 2385 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2386 aa = a->a; 2387 ai = a->i; 2388 aj = a->j; 2389 2390 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2391 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2392 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2393 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2394 for (i=0; i<m; i++) { 2395 ncols = ai[1] - ai[0]; ai++; 2396 if (ncols == A->cmap->n) { /* row is dense */ 2397 x[i] = *aa; if (idx) idx[i] = 0; 2398 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2399 x[i] = 0.0; 2400 if (idx) { 2401 idx[i] = 0; /* in case ncols is zero */ 2402 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2403 if (aj[j] > j) { 2404 idx[i] = j; 2405 break; 2406 } 2407 } 2408 } 2409 } 2410 for (j=0; j<ncols; j++){ 2411 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2412 aa++; aj++; 2413 } 2414 } 2415 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2416 PetscFunctionReturn(0); 2417 } 2418 2419 #undef __FUNCT__ 2420 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2421 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2422 { 2423 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2424 PetscErrorCode ierr; 2425 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2426 PetscReal atmp; 2427 PetscScalar *x; 2428 MatScalar *aa; 2429 2430 PetscFunctionBegin; 2431 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2432 aa = a->a; 2433 ai = a->i; 2434 aj = a->j; 2435 2436 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2437 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2438 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2439 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2440 for (i=0; i<m; i++) { 2441 ncols = ai[1] - ai[0]; ai++; 2442 if (ncols) { 2443 /* Get first nonzero */ 2444 for(j = 0; j < ncols; j++) { 2445 atmp = PetscAbsScalar(aa[j]); 2446 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2447 } 2448 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2449 } else { 2450 x[i] = 0.0; if (idx) idx[i] = 0; 2451 } 2452 for(j = 0; j < ncols; j++) { 2453 atmp = PetscAbsScalar(*aa); 2454 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2455 aa++; aj++; 2456 } 2457 } 2458 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2459 PetscFunctionReturn(0); 2460 } 2461 2462 #undef __FUNCT__ 2463 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2464 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2465 { 2466 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2467 PetscErrorCode ierr; 2468 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2469 PetscScalar *x; 2470 MatScalar *aa; 2471 2472 PetscFunctionBegin; 2473 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2474 aa = a->a; 2475 ai = a->i; 2476 aj = a->j; 2477 2478 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2479 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2480 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2481 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2482 for (i=0; i<m; i++) { 2483 ncols = ai[1] - ai[0]; ai++; 2484 if (ncols == A->cmap->n) { /* row is dense */ 2485 x[i] = *aa; if (idx) idx[i] = 0; 2486 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2487 x[i] = 0.0; 2488 if (idx) { /* find first implicit 0.0 in the row */ 2489 idx[i] = 0; /* in case ncols is zero */ 2490 for (j=0;j<ncols;j++) { 2491 if (aj[j] > j) { 2492 idx[i] = j; 2493 break; 2494 } 2495 } 2496 } 2497 } 2498 for (j=0; j<ncols; j++){ 2499 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2500 aa++; aj++; 2501 } 2502 } 2503 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2504 PetscFunctionReturn(0); 2505 } 2506 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2507 /* -------------------------------------------------------------------*/ 2508 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2509 MatGetRow_SeqAIJ, 2510 MatRestoreRow_SeqAIJ, 2511 MatMult_SeqAIJ, 2512 /* 4*/ MatMultAdd_SeqAIJ, 2513 MatMultTranspose_SeqAIJ, 2514 MatMultTransposeAdd_SeqAIJ, 2515 0, 2516 0, 2517 0, 2518 /*10*/ 0, 2519 MatLUFactor_SeqAIJ, 2520 0, 2521 MatSOR_SeqAIJ, 2522 MatTranspose_SeqAIJ, 2523 /*15*/ MatGetInfo_SeqAIJ, 2524 MatEqual_SeqAIJ, 2525 MatGetDiagonal_SeqAIJ, 2526 MatDiagonalScale_SeqAIJ, 2527 MatNorm_SeqAIJ, 2528 /*20*/ 0, 2529 MatAssemblyEnd_SeqAIJ, 2530 MatSetOption_SeqAIJ, 2531 MatZeroEntries_SeqAIJ, 2532 /*24*/ MatZeroRows_SeqAIJ, 2533 0, 2534 0, 2535 0, 2536 0, 2537 /*29*/ MatSetUpPreallocation_SeqAIJ, 2538 0, 2539 0, 2540 MatGetArray_SeqAIJ, 2541 MatRestoreArray_SeqAIJ, 2542 /*34*/ MatDuplicate_SeqAIJ, 2543 0, 2544 0, 2545 MatILUFactor_SeqAIJ, 2546 0, 2547 /*39*/ MatAXPY_SeqAIJ, 2548 MatGetSubMatrices_SeqAIJ, 2549 MatIncreaseOverlap_SeqAIJ, 2550 MatGetValues_SeqAIJ, 2551 MatCopy_SeqAIJ, 2552 /*44*/ MatGetRowMax_SeqAIJ, 2553 MatScale_SeqAIJ, 2554 0, 2555 MatDiagonalSet_SeqAIJ, 2556 0, 2557 /*49*/ MatSetBlockSize_SeqAIJ, 2558 MatGetRowIJ_SeqAIJ, 2559 MatRestoreRowIJ_SeqAIJ, 2560 MatGetColumnIJ_SeqAIJ, 2561 MatRestoreColumnIJ_SeqAIJ, 2562 /*54*/ MatFDColoringCreate_SeqAIJ, 2563 0, 2564 0, 2565 MatPermute_SeqAIJ, 2566 0, 2567 /*59*/ 0, 2568 MatDestroy_SeqAIJ, 2569 MatView_SeqAIJ, 2570 0, 2571 0, 2572 /*64*/ 0, 2573 0, 2574 0, 2575 0, 2576 0, 2577 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2578 MatGetRowMinAbs_SeqAIJ, 2579 0, 2580 MatSetColoring_SeqAIJ, 2581 #if defined(PETSC_HAVE_ADIC) 2582 MatSetValuesAdic_SeqAIJ, 2583 #else 2584 0, 2585 #endif 2586 /*74*/ MatSetValuesAdifor_SeqAIJ, 2587 MatFDColoringApply_AIJ, 2588 0, 2589 0, 2590 0, 2591 /*79*/ 0, 2592 0, 2593 0, 2594 0, 2595 MatLoad_SeqAIJ, 2596 /*84*/ MatIsSymmetric_SeqAIJ, 2597 MatIsHermitian_SeqAIJ, 2598 0, 2599 0, 2600 0, 2601 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2602 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2603 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2604 MatPtAP_Basic, 2605 MatPtAPSymbolic_SeqAIJ, 2606 /*94*/ MatPtAPNumeric_SeqAIJ, 2607 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2608 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2609 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2610 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2611 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2612 0, 2613 0, 2614 MatConjugate_SeqAIJ, 2615 0, 2616 /*104*/MatSetValuesRow_SeqAIJ, 2617 MatRealPart_SeqAIJ, 2618 MatImaginaryPart_SeqAIJ, 2619 0, 2620 0, 2621 /*109*/0, 2622 0, 2623 MatGetRowMin_SeqAIJ, 2624 0, 2625 MatMissingDiagonal_SeqAIJ, 2626 /*114*/0, 2627 0, 2628 0, 2629 0, 2630 0, 2631 /*119*/0, 2632 0, 2633 0, 2634 0, 2635 MatLoadnew_SeqAIJ 2636 }; 2637 2638 EXTERN_C_BEGIN 2639 #undef __FUNCT__ 2640 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2641 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2642 { 2643 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2644 PetscInt i,nz,n; 2645 2646 PetscFunctionBegin; 2647 2648 nz = aij->maxnz; 2649 n = mat->rmap->n; 2650 for (i=0; i<nz; i++) { 2651 aij->j[i] = indices[i]; 2652 } 2653 aij->nz = nz; 2654 for (i=0; i<n; i++) { 2655 aij->ilen[i] = aij->imax[i]; 2656 } 2657 2658 PetscFunctionReturn(0); 2659 } 2660 EXTERN_C_END 2661 2662 #undef __FUNCT__ 2663 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2664 /*@ 2665 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2666 in the matrix. 2667 2668 Input Parameters: 2669 + mat - the SeqAIJ matrix 2670 - indices - the column indices 2671 2672 Level: advanced 2673 2674 Notes: 2675 This can be called if you have precomputed the nonzero structure of the 2676 matrix and want to provide it to the matrix object to improve the performance 2677 of the MatSetValues() operation. 2678 2679 You MUST have set the correct numbers of nonzeros per row in the call to 2680 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2681 2682 MUST be called before any calls to MatSetValues(); 2683 2684 The indices should start with zero, not one. 2685 2686 @*/ 2687 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2688 { 2689 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2690 2691 PetscFunctionBegin; 2692 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2693 PetscValidPointer(indices,2); 2694 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2695 if (f) { 2696 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2697 } else { 2698 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2699 } 2700 PetscFunctionReturn(0); 2701 } 2702 2703 /* ----------------------------------------------------------------------------------------*/ 2704 2705 EXTERN_C_BEGIN 2706 #undef __FUNCT__ 2707 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2708 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2709 { 2710 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2711 PetscErrorCode ierr; 2712 size_t nz = aij->i[mat->rmap->n]; 2713 2714 PetscFunctionBegin; 2715 if (aij->nonew != 1) { 2716 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2717 } 2718 2719 /* allocate space for values if not already there */ 2720 if (!aij->saved_values) { 2721 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2722 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2723 } 2724 2725 /* copy values over */ 2726 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2727 PetscFunctionReturn(0); 2728 } 2729 EXTERN_C_END 2730 2731 #undef __FUNCT__ 2732 #define __FUNCT__ "MatStoreValues" 2733 /*@ 2734 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2735 example, reuse of the linear part of a Jacobian, while recomputing the 2736 nonlinear portion. 2737 2738 Collect on Mat 2739 2740 Input Parameters: 2741 . mat - the matrix (currently only AIJ matrices support this option) 2742 2743 Level: advanced 2744 2745 Common Usage, with SNESSolve(): 2746 $ Create Jacobian matrix 2747 $ Set linear terms into matrix 2748 $ Apply boundary conditions to matrix, at this time matrix must have 2749 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2750 $ boundary conditions again will not change the nonzero structure 2751 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2752 $ ierr = MatStoreValues(mat); 2753 $ Call SNESSetJacobian() with matrix 2754 $ In your Jacobian routine 2755 $ ierr = MatRetrieveValues(mat); 2756 $ Set nonlinear terms in matrix 2757 2758 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2759 $ // build linear portion of Jacobian 2760 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2761 $ ierr = MatStoreValues(mat); 2762 $ loop over nonlinear iterations 2763 $ ierr = MatRetrieveValues(mat); 2764 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2765 $ // call MatAssemblyBegin/End() on matrix 2766 $ Solve linear system with Jacobian 2767 $ endloop 2768 2769 Notes: 2770 Matrix must already be assemblied before calling this routine 2771 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2772 calling this routine. 2773 2774 When this is called multiple times it overwrites the previous set of stored values 2775 and does not allocated additional space. 2776 2777 .seealso: MatRetrieveValues() 2778 2779 @*/ 2780 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2781 { 2782 PetscErrorCode ierr,(*f)(Mat); 2783 2784 PetscFunctionBegin; 2785 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2786 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2787 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2788 2789 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2790 if (f) { 2791 ierr = (*f)(mat);CHKERRQ(ierr); 2792 } else { 2793 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2794 } 2795 PetscFunctionReturn(0); 2796 } 2797 2798 EXTERN_C_BEGIN 2799 #undef __FUNCT__ 2800 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2801 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2802 { 2803 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2804 PetscErrorCode ierr; 2805 PetscInt nz = aij->i[mat->rmap->n]; 2806 2807 PetscFunctionBegin; 2808 if (aij->nonew != 1) { 2809 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2810 } 2811 if (!aij->saved_values) { 2812 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2813 } 2814 /* copy values over */ 2815 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2816 PetscFunctionReturn(0); 2817 } 2818 EXTERN_C_END 2819 2820 #undef __FUNCT__ 2821 #define __FUNCT__ "MatRetrieveValues" 2822 /*@ 2823 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2824 example, reuse of the linear part of a Jacobian, while recomputing the 2825 nonlinear portion. 2826 2827 Collect on Mat 2828 2829 Input Parameters: 2830 . mat - the matrix (currently on AIJ matrices support this option) 2831 2832 Level: advanced 2833 2834 .seealso: MatStoreValues() 2835 2836 @*/ 2837 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2838 { 2839 PetscErrorCode ierr,(*f)(Mat); 2840 2841 PetscFunctionBegin; 2842 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2843 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2844 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2845 2846 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2847 if (f) { 2848 ierr = (*f)(mat);CHKERRQ(ierr); 2849 } else { 2850 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2851 } 2852 PetscFunctionReturn(0); 2853 } 2854 2855 2856 /* --------------------------------------------------------------------------------*/ 2857 #undef __FUNCT__ 2858 #define __FUNCT__ "MatCreateSeqAIJ" 2859 /*@C 2860 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2861 (the default parallel PETSc format). For good matrix assembly performance 2862 the user should preallocate the matrix storage by setting the parameter nz 2863 (or the array nnz). By setting these parameters accurately, performance 2864 during matrix assembly can be increased by more than a factor of 50. 2865 2866 Collective on MPI_Comm 2867 2868 Input Parameters: 2869 + comm - MPI communicator, set to PETSC_COMM_SELF 2870 . m - number of rows 2871 . n - number of columns 2872 . nz - number of nonzeros per row (same for all rows) 2873 - nnz - array containing the number of nonzeros in the various rows 2874 (possibly different for each row) or PETSC_NULL 2875 2876 Output Parameter: 2877 . A - the matrix 2878 2879 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2880 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2881 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2882 2883 Notes: 2884 If nnz is given then nz is ignored 2885 2886 The AIJ format (also called the Yale sparse matrix format or 2887 compressed row storage), is fully compatible with standard Fortran 77 2888 storage. That is, the stored row and column indices can begin at 2889 either one (as in Fortran) or zero. See the users' manual for details. 2890 2891 Specify the preallocated storage with either nz or nnz (not both). 2892 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2893 allocation. For large problems you MUST preallocate memory or you 2894 will get TERRIBLE performance, see the users' manual chapter on matrices. 2895 2896 By default, this format uses inodes (identical nodes) when possible, to 2897 improve numerical efficiency of matrix-vector products and solves. We 2898 search for consecutive rows with the same nonzero structure, thereby 2899 reusing matrix information to achieve increased efficiency. 2900 2901 Options Database Keys: 2902 + -mat_no_inode - Do not use inodes 2903 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2904 - -mat_aij_oneindex - Internally use indexing starting at 1 2905 rather than 0. Note that when calling MatSetValues(), 2906 the user still MUST index entries starting at 0! 2907 2908 Level: intermediate 2909 2910 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2911 2912 @*/ 2913 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2914 { 2915 PetscErrorCode ierr; 2916 2917 PetscFunctionBegin; 2918 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2919 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2920 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2921 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2922 PetscFunctionReturn(0); 2923 } 2924 2925 #undef __FUNCT__ 2926 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2927 /*@C 2928 MatSeqAIJSetPreallocation - For good matrix assembly performance 2929 the user should preallocate the matrix storage by setting the parameter nz 2930 (or the array nnz). By setting these parameters accurately, performance 2931 during matrix assembly can be increased by more than a factor of 50. 2932 2933 Collective on MPI_Comm 2934 2935 Input Parameters: 2936 + B - The matrix-free 2937 . nz - number of nonzeros per row (same for all rows) 2938 - nnz - array containing the number of nonzeros in the various rows 2939 (possibly different for each row) or PETSC_NULL 2940 2941 Notes: 2942 If nnz is given then nz is ignored 2943 2944 The AIJ format (also called the Yale sparse matrix format or 2945 compressed row storage), is fully compatible with standard Fortran 77 2946 storage. That is, the stored row and column indices can begin at 2947 either one (as in Fortran) or zero. See the users' manual for details. 2948 2949 Specify the preallocated storage with either nz or nnz (not both). 2950 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2951 allocation. For large problems you MUST preallocate memory or you 2952 will get TERRIBLE performance, see the users' manual chapter on matrices. 2953 2954 You can call MatGetInfo() to get information on how effective the preallocation was; 2955 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2956 You can also run with the option -info and look for messages with the string 2957 malloc in them to see if additional memory allocation was needed. 2958 2959 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2960 entries or columns indices 2961 2962 By default, this format uses inodes (identical nodes) when possible, to 2963 improve numerical efficiency of matrix-vector products and solves. We 2964 search for consecutive rows with the same nonzero structure, thereby 2965 reusing matrix information to achieve increased efficiency. 2966 2967 Options Database Keys: 2968 + -mat_no_inode - Do not use inodes 2969 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2970 - -mat_aij_oneindex - Internally use indexing starting at 1 2971 rather than 0. Note that when calling MatSetValues(), 2972 the user still MUST index entries starting at 0! 2973 2974 Level: intermediate 2975 2976 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2977 2978 @*/ 2979 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2980 { 2981 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2982 2983 PetscFunctionBegin; 2984 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2985 if (f) { 2986 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2987 } 2988 PetscFunctionReturn(0); 2989 } 2990 2991 EXTERN_C_BEGIN 2992 #undef __FUNCT__ 2993 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2994 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2995 { 2996 Mat_SeqAIJ *b; 2997 PetscTruth skipallocation = PETSC_FALSE; 2998 PetscErrorCode ierr; 2999 PetscInt i; 3000 3001 PetscFunctionBegin; 3002 3003 if (nz == MAT_SKIP_ALLOCATION) { 3004 skipallocation = PETSC_TRUE; 3005 nz = 0; 3006 } 3007 3008 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3009 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3010 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3011 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3012 3013 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3014 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3015 if (nnz) { 3016 for (i=0; i<B->rmap->n; i++) { 3017 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]); 3018 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); 3019 } 3020 } 3021 3022 B->preallocated = PETSC_TRUE; 3023 b = (Mat_SeqAIJ*)B->data; 3024 3025 if (!skipallocation) { 3026 if (!b->imax) { 3027 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3028 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3029 } 3030 if (!nnz) { 3031 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3032 else if (nz <= 0) nz = 1; 3033 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3034 nz = nz*B->rmap->n; 3035 } else { 3036 nz = 0; 3037 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3038 } 3039 /* b->ilen will count nonzeros in each row so far. */ 3040 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3041 3042 /* allocate the matrix space */ 3043 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3044 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3045 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3046 b->i[0] = 0; 3047 for (i=1; i<B->rmap->n+1; i++) { 3048 b->i[i] = b->i[i-1] + b->imax[i-1]; 3049 } 3050 b->singlemalloc = PETSC_TRUE; 3051 b->free_a = PETSC_TRUE; 3052 b->free_ij = PETSC_TRUE; 3053 } else { 3054 b->free_a = PETSC_FALSE; 3055 b->free_ij = PETSC_FALSE; 3056 } 3057 3058 b->nz = 0; 3059 b->maxnz = nz; 3060 B->info.nz_unneeded = (double)b->maxnz; 3061 PetscFunctionReturn(0); 3062 } 3063 EXTERN_C_END 3064 3065 #undef __FUNCT__ 3066 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3067 /*@ 3068 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3069 3070 Input Parameters: 3071 + B - the matrix 3072 . i - the indices into j for the start of each row (starts with zero) 3073 . j - the column indices for each row (starts with zero) these must be sorted for each row 3074 - v - optional values in the matrix 3075 3076 Level: developer 3077 3078 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3079 3080 .keywords: matrix, aij, compressed row, sparse, sequential 3081 3082 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3083 @*/ 3084 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3085 { 3086 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3087 PetscErrorCode ierr; 3088 3089 PetscFunctionBegin; 3090 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3091 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3092 if (f) { 3093 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3094 } 3095 PetscFunctionReturn(0); 3096 } 3097 3098 EXTERN_C_BEGIN 3099 #undef __FUNCT__ 3100 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3101 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3102 { 3103 PetscInt i; 3104 PetscInt m,n; 3105 PetscInt nz; 3106 PetscInt *nnz, nz_max = 0; 3107 PetscScalar *values; 3108 PetscErrorCode ierr; 3109 3110 PetscFunctionBegin; 3111 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3112 3113 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3114 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3115 for(i = 0; i < m; i++) { 3116 nz = Ii[i+1]- Ii[i]; 3117 nz_max = PetscMax(nz_max, nz); 3118 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3119 nnz[i] = nz; 3120 } 3121 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3122 ierr = PetscFree(nnz);CHKERRQ(ierr); 3123 3124 if (v) { 3125 values = (PetscScalar*) v; 3126 } else { 3127 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3128 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3129 } 3130 3131 for(i = 0; i < m; i++) { 3132 nz = Ii[i+1] - Ii[i]; 3133 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3134 } 3135 3136 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3137 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3138 3139 if (!v) { 3140 ierr = PetscFree(values);CHKERRQ(ierr); 3141 } 3142 PetscFunctionReturn(0); 3143 } 3144 EXTERN_C_END 3145 3146 #include "../src/mat/impls/dense/seq/dense.h" 3147 #include "private/petscaxpy.h" 3148 3149 #undef __FUNCT__ 3150 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3151 /* 3152 Computes (B'*A')' since computing B*A directly is untenable 3153 3154 n p p 3155 ( ) ( ) ( ) 3156 m ( A ) * n ( B ) = m ( C ) 3157 ( ) ( ) ( ) 3158 3159 */ 3160 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3161 { 3162 PetscErrorCode ierr; 3163 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3164 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3165 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3166 PetscInt i,n,m,q,p; 3167 const PetscInt *ii,*idx; 3168 const PetscScalar *b,*a,*a_q; 3169 PetscScalar *c,*c_q; 3170 3171 PetscFunctionBegin; 3172 m = A->rmap->n; 3173 n = A->cmap->n; 3174 p = B->cmap->n; 3175 a = sub_a->v; 3176 b = sub_b->a; 3177 c = sub_c->v; 3178 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3179 3180 ii = sub_b->i; 3181 idx = sub_b->j; 3182 for (i=0; i<n; i++) { 3183 q = ii[i+1] - ii[i]; 3184 while (q-->0) { 3185 c_q = c + m*(*idx); 3186 a_q = a + m*i; 3187 PetscAXPY(c_q,*b,a_q,m); 3188 idx++; 3189 b++; 3190 } 3191 } 3192 PetscFunctionReturn(0); 3193 } 3194 3195 #undef __FUNCT__ 3196 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3197 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3198 { 3199 PetscErrorCode ierr; 3200 PetscInt m=A->rmap->n,n=B->cmap->n; 3201 Mat Cmat; 3202 3203 PetscFunctionBegin; 3204 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); 3205 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3206 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3207 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3208 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3209 Cmat->assembled = PETSC_TRUE; 3210 *C = Cmat; 3211 PetscFunctionReturn(0); 3212 } 3213 3214 /* ----------------------------------------------------------------*/ 3215 #undef __FUNCT__ 3216 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3217 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3218 { 3219 PetscErrorCode ierr; 3220 3221 PetscFunctionBegin; 3222 if (scall == MAT_INITIAL_MATRIX){ 3223 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3224 } 3225 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3226 PetscFunctionReturn(0); 3227 } 3228 3229 3230 /*MC 3231 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3232 based on compressed sparse row format. 3233 3234 Options Database Keys: 3235 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3236 3237 Level: beginner 3238 3239 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3240 M*/ 3241 3242 EXTERN_C_BEGIN 3243 #if defined(PETSC_HAVE_PASTIX) 3244 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3245 #endif 3246 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3247 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3248 #endif 3249 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3250 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3251 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3252 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3253 #if defined(PETSC_HAVE_MUMPS) 3254 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3255 #endif 3256 #if defined(PETSC_HAVE_SUPERLU) 3257 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3258 #endif 3259 #if defined(PETSC_HAVE_SUPERLU_DIST) 3260 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3261 #endif 3262 #if defined(PETSC_HAVE_SPOOLES) 3263 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3264 #endif 3265 #if defined(PETSC_HAVE_UMFPACK) 3266 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3267 #endif 3268 #if defined(PETSC_HAVE_CHOLMOD) 3269 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3270 #endif 3271 #if defined(PETSC_HAVE_LUSOL) 3272 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3273 #endif 3274 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3275 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3276 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3277 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3278 #endif 3279 EXTERN_C_END 3280 3281 3282 EXTERN_C_BEGIN 3283 #undef __FUNCT__ 3284 #define __FUNCT__ "MatCreate_SeqAIJ" 3285 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3286 { 3287 Mat_SeqAIJ *b; 3288 PetscErrorCode ierr; 3289 PetscMPIInt size; 3290 3291 PetscFunctionBegin; 3292 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3293 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3294 3295 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3296 B->data = (void*)b; 3297 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3298 B->mapping = 0; 3299 b->row = 0; 3300 b->col = 0; 3301 b->icol = 0; 3302 b->reallocs = 0; 3303 b->ignorezeroentries = PETSC_FALSE; 3304 b->roworiented = PETSC_TRUE; 3305 b->nonew = 0; 3306 b->diag = 0; 3307 b->solve_work = 0; 3308 B->spptr = 0; 3309 b->saved_values = 0; 3310 b->idiag = 0; 3311 b->mdiag = 0; 3312 b->ssor_work = 0; 3313 b->omega = 1.0; 3314 b->fshift = 0.0; 3315 b->idiagvalid = PETSC_FALSE; 3316 b->keepnonzeropattern = PETSC_FALSE; 3317 b->xtoy = 0; 3318 b->XtoY = 0; 3319 b->compressedrow.use = PETSC_FALSE; 3320 b->compressedrow.nrows = B->rmap->n; 3321 b->compressedrow.i = PETSC_NULL; 3322 b->compressedrow.rindex = PETSC_NULL; 3323 b->compressedrow.checked = PETSC_FALSE; 3324 B->same_nonzero = PETSC_FALSE; 3325 3326 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3327 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3328 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3329 "MatGetFactor_seqaij_matlab", 3330 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3331 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3332 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3333 #endif 3334 #if defined(PETSC_HAVE_PASTIX) 3335 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3336 "MatGetFactor_seqaij_pastix", 3337 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3338 #endif 3339 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3340 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3341 "MatGetFactor_seqaij_essl", 3342 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3343 #endif 3344 #if defined(PETSC_HAVE_SUPERLU) 3345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3346 "MatGetFactor_seqaij_superlu", 3347 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3348 #endif 3349 #if defined(PETSC_HAVE_SUPERLU_DIST) 3350 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3351 "MatGetFactor_seqaij_superlu_dist", 3352 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3353 #endif 3354 #if defined(PETSC_HAVE_SPOOLES) 3355 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3356 "MatGetFactor_seqaij_spooles", 3357 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3358 #endif 3359 #if defined(PETSC_HAVE_MUMPS) 3360 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3361 "MatGetFactor_aij_mumps", 3362 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3363 #endif 3364 #if defined(PETSC_HAVE_UMFPACK) 3365 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3366 "MatGetFactor_seqaij_umfpack", 3367 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3368 #endif 3369 #if defined(PETSC_HAVE_CHOLMOD) 3370 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3371 "MatGetFactor_seqaij_cholmod", 3372 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3373 #endif 3374 #if defined(PETSC_HAVE_LUSOL) 3375 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3376 "MatGetFactor_seqaij_lusol", 3377 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3378 #endif 3379 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3380 "MatGetFactor_seqaij_petsc", 3381 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3382 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3383 "MatGetFactorAvailable_seqaij_petsc", 3384 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3385 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3386 "MatGetFactor_seqaij_bas", 3387 MatGetFactor_seqaij_bas); 3388 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3389 "MatSeqAIJSetColumnIndices_SeqAIJ", 3390 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3391 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3392 "MatStoreValues_SeqAIJ", 3393 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3394 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3395 "MatRetrieveValues_SeqAIJ", 3396 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3397 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3398 "MatConvert_SeqAIJ_SeqSBAIJ", 3399 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3400 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3401 "MatConvert_SeqAIJ_SeqBAIJ", 3402 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3403 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3404 "MatConvert_SeqAIJ_SeqCSRPERM", 3405 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3406 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3407 "MatConvert_SeqAIJ_SeqCRL", 3408 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3409 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3410 "MatIsTranspose_SeqAIJ", 3411 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3412 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3413 "MatIsHermitianTranspose_SeqAIJ", 3414 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3415 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3416 "MatSeqAIJSetPreallocation_SeqAIJ", 3417 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3418 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3419 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3420 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3421 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3422 "MatReorderForNonzeroDiagonal_SeqAIJ", 3423 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3424 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3425 "MatMatMult_SeqDense_SeqAIJ", 3426 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3427 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3428 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3429 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3430 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3431 "MatMatMultNumeric_SeqDense_SeqAIJ", 3432 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3433 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3434 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3435 PetscFunctionReturn(0); 3436 } 3437 EXTERN_C_END 3438 3439 #undef __FUNCT__ 3440 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3441 /* 3442 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3443 */ 3444 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3445 { 3446 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3447 PetscErrorCode ierr; 3448 PetscInt i,m = A->rmap->n; 3449 3450 PetscFunctionBegin; 3451 c = (Mat_SeqAIJ*)C->data; 3452 3453 C->factortype = A->factortype; 3454 c->row = 0; 3455 c->col = 0; 3456 c->icol = 0; 3457 c->reallocs = 0; 3458 3459 C->assembled = PETSC_TRUE; 3460 3461 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3462 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3463 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3464 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3465 3466 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3467 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3468 for (i=0; i<m; i++) { 3469 c->imax[i] = a->imax[i]; 3470 c->ilen[i] = a->ilen[i]; 3471 } 3472 3473 /* allocate the matrix space */ 3474 if (mallocmatspace){ 3475 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3476 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3477 c->singlemalloc = PETSC_TRUE; 3478 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3479 if (m > 0) { 3480 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3481 if (cpvalues == MAT_COPY_VALUES) { 3482 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3483 } else { 3484 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3485 } 3486 } 3487 } 3488 3489 c->ignorezeroentries = a->ignorezeroentries; 3490 c->roworiented = a->roworiented; 3491 c->nonew = a->nonew; 3492 if (a->diag) { 3493 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3494 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3495 for (i=0; i<m; i++) { 3496 c->diag[i] = a->diag[i]; 3497 } 3498 } else c->diag = 0; 3499 c->solve_work = 0; 3500 c->saved_values = 0; 3501 c->idiag = 0; 3502 c->ssor_work = 0; 3503 c->keepnonzeropattern = a->keepnonzeropattern; 3504 c->free_a = PETSC_TRUE; 3505 c->free_ij = PETSC_TRUE; 3506 c->xtoy = 0; 3507 c->XtoY = 0; 3508 3509 c->nz = a->nz; 3510 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3511 C->preallocated = PETSC_TRUE; 3512 3513 c->compressedrow.use = a->compressedrow.use; 3514 c->compressedrow.nrows = a->compressedrow.nrows; 3515 c->compressedrow.checked = a->compressedrow.checked; 3516 if (a->compressedrow.checked && a->compressedrow.use){ 3517 i = a->compressedrow.nrows; 3518 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3519 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3520 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3521 } else { 3522 c->compressedrow.use = PETSC_FALSE; 3523 c->compressedrow.i = PETSC_NULL; 3524 c->compressedrow.rindex = PETSC_NULL; 3525 } 3526 C->same_nonzero = A->same_nonzero; 3527 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3528 3529 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3530 PetscFunctionReturn(0); 3531 } 3532 3533 #undef __FUNCT__ 3534 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3535 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3536 { 3537 PetscErrorCode ierr; 3538 3539 PetscFunctionBegin; 3540 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3541 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3542 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3543 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3544 PetscFunctionReturn(0); 3545 } 3546 3547 #undef __FUNCT__ 3548 #define __FUNCT__ "MatLoad_SeqAIJ" 3549 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, const MatType type,Mat *A) 3550 { 3551 Mat_SeqAIJ *a; 3552 Mat B; 3553 PetscErrorCode ierr; 3554 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3555 int fd; 3556 PetscMPIInt size; 3557 MPI_Comm comm; 3558 3559 PetscFunctionBegin; 3560 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3561 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3562 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3563 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3564 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3565 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3566 M = header[1]; N = header[2]; nz = header[3]; 3567 3568 if (nz < 0) { 3569 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3570 } 3571 3572 /* read in row lengths */ 3573 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3574 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3575 3576 /* check if sum of rowlengths is same as nz */ 3577 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3578 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); 3579 3580 /* create our matrix */ 3581 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3582 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3583 ierr = MatSetType(B,type);CHKERRQ(ierr); 3584 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3585 a = (Mat_SeqAIJ*)B->data; 3586 3587 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3588 3589 /* read in nonzero values */ 3590 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3591 3592 /* set matrix "i" values */ 3593 a->i[0] = 0; 3594 for (i=1; i<= M; i++) { 3595 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3596 a->ilen[i-1] = rowlengths[i-1]; 3597 } 3598 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3599 3600 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3601 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3602 *A = B; 3603 PetscFunctionReturn(0); 3604 } 3605 3606 #undef __FUNCT__ 3607 #define __FUNCT__ "MatLoadnew_SeqAIJ" 3608 PetscErrorCode MatLoadnew_SeqAIJ(PetscViewer viewer,Mat newMat) 3609 { 3610 Mat_SeqAIJ *a; 3611 PetscErrorCode ierr; 3612 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3613 int fd; 3614 PetscMPIInt size; 3615 MPI_Comm comm; 3616 3617 PetscFunctionBegin; 3618 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3619 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3620 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3621 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3622 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3623 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3624 M = header[1]; N = header[2]; nz = header[3]; 3625 3626 if (nz < 0) { 3627 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3628 } 3629 3630 /* read in row lengths */ 3631 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3632 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3633 3634 /* check if sum of rowlengths is same as nz */ 3635 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3636 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); 3637 3638 /* set global size if not set already*/ 3639 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3640 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3641 } else { 3642 /* if sizes and type are already set, check if the vector global sizes are correct */ 3643 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3644 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); 3645 } 3646 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3647 a = (Mat_SeqAIJ*)newMat->data; 3648 3649 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3650 3651 /* read in nonzero values */ 3652 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3653 3654 /* set matrix "i" values */ 3655 a->i[0] = 0; 3656 for (i=1; i<= M; i++) { 3657 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3658 a->ilen[i-1] = rowlengths[i-1]; 3659 } 3660 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3661 3662 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3663 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3664 PetscFunctionReturn(0); 3665 } 3666 3667 #undef __FUNCT__ 3668 #define __FUNCT__ "MatEqual_SeqAIJ" 3669 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3670 { 3671 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3672 PetscErrorCode ierr; 3673 #if defined(PETSC_USE_COMPLEX) 3674 PetscInt k; 3675 #endif 3676 3677 PetscFunctionBegin; 3678 /* If the matrix dimensions are not equal,or no of nonzeros */ 3679 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3680 *flg = PETSC_FALSE; 3681 PetscFunctionReturn(0); 3682 } 3683 3684 /* if the a->i are the same */ 3685 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3686 if (!*flg) PetscFunctionReturn(0); 3687 3688 /* if a->j are the same */ 3689 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3690 if (!*flg) PetscFunctionReturn(0); 3691 3692 /* if a->a are the same */ 3693 #if defined(PETSC_USE_COMPLEX) 3694 for (k=0; k<a->nz; k++){ 3695 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3696 *flg = PETSC_FALSE; 3697 PetscFunctionReturn(0); 3698 } 3699 } 3700 #else 3701 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3702 #endif 3703 PetscFunctionReturn(0); 3704 } 3705 3706 #undef __FUNCT__ 3707 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3708 /*@ 3709 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3710 provided by the user. 3711 3712 Collective on MPI_Comm 3713 3714 Input Parameters: 3715 + comm - must be an MPI communicator of size 1 3716 . m - number of rows 3717 . n - number of columns 3718 . i - row indices 3719 . j - column indices 3720 - a - matrix values 3721 3722 Output Parameter: 3723 . mat - the matrix 3724 3725 Level: intermediate 3726 3727 Notes: 3728 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3729 once the matrix is destroyed 3730 3731 You cannot set new nonzero locations into this matrix, that will generate an error. 3732 3733 The i and j indices are 0 based 3734 3735 The format which is used for the sparse matrix input, is equivalent to a 3736 row-major ordering.. i.e for the following matrix, the input data expected is 3737 as shown: 3738 3739 1 0 0 3740 2 0 3 3741 4 5 6 3742 3743 i = {0,1,3,6} [size = nrow+1 = 3+1] 3744 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3745 v = {1,2,3,4,5,6} [size = nz = 6] 3746 3747 3748 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3749 3750 @*/ 3751 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3752 { 3753 PetscErrorCode ierr; 3754 PetscInt ii; 3755 Mat_SeqAIJ *aij; 3756 #if defined(PETSC_USE_DEBUG) 3757 PetscInt jj; 3758 #endif 3759 3760 PetscFunctionBegin; 3761 if (i[0]) { 3762 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3763 } 3764 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3765 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3766 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3767 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3768 aij = (Mat_SeqAIJ*)(*mat)->data; 3769 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3770 3771 aij->i = i; 3772 aij->j = j; 3773 aij->a = a; 3774 aij->singlemalloc = PETSC_FALSE; 3775 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3776 aij->free_a = PETSC_FALSE; 3777 aij->free_ij = PETSC_FALSE; 3778 3779 for (ii=0; ii<m; ii++) { 3780 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3781 #if defined(PETSC_USE_DEBUG) 3782 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]); 3783 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3784 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); 3785 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); 3786 } 3787 #endif 3788 } 3789 #if defined(PETSC_USE_DEBUG) 3790 for (ii=0; ii<aij->i[m]; ii++) { 3791 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3792 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]); 3793 } 3794 #endif 3795 3796 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3797 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3798 PetscFunctionReturn(0); 3799 } 3800 3801 #undef __FUNCT__ 3802 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3803 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3804 { 3805 PetscErrorCode ierr; 3806 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3807 3808 PetscFunctionBegin; 3809 if (coloring->ctype == IS_COLORING_GLOBAL) { 3810 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3811 a->coloring = coloring; 3812 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3813 PetscInt i,*larray; 3814 ISColoring ocoloring; 3815 ISColoringValue *colors; 3816 3817 /* set coloring for diagonal portion */ 3818 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3819 for (i=0; i<A->cmap->n; i++) { 3820 larray[i] = i; 3821 } 3822 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3823 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3824 for (i=0; i<A->cmap->n; i++) { 3825 colors[i] = coloring->colors[larray[i]]; 3826 } 3827 ierr = PetscFree(larray);CHKERRQ(ierr); 3828 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3829 a->coloring = ocoloring; 3830 } 3831 PetscFunctionReturn(0); 3832 } 3833 3834 #if defined(PETSC_HAVE_ADIC) 3835 EXTERN_C_BEGIN 3836 #include "adic/ad_utils.h" 3837 EXTERN_C_END 3838 3839 #undef __FUNCT__ 3840 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3841 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3842 { 3843 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3844 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3845 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3846 ISColoringValue *color; 3847 3848 PetscFunctionBegin; 3849 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3850 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3851 color = a->coloring->colors; 3852 /* loop over rows */ 3853 for (i=0; i<m; i++) { 3854 nz = ii[i+1] - ii[i]; 3855 /* loop over columns putting computed value into matrix */ 3856 for (j=0; j<nz; j++) { 3857 *v++ = values[color[*jj++]]; 3858 } 3859 values += nlen; /* jump to next row of derivatives */ 3860 } 3861 PetscFunctionReturn(0); 3862 } 3863 #endif 3864 3865 #undef __FUNCT__ 3866 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3867 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3868 { 3869 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3870 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3871 MatScalar *v = a->a; 3872 PetscScalar *values = (PetscScalar *)advalues; 3873 ISColoringValue *color; 3874 3875 PetscFunctionBegin; 3876 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3877 color = a->coloring->colors; 3878 /* loop over rows */ 3879 for (i=0; i<m; i++) { 3880 nz = ii[i+1] - ii[i]; 3881 /* loop over columns putting computed value into matrix */ 3882 for (j=0; j<nz; j++) { 3883 *v++ = values[color[*jj++]]; 3884 } 3885 values += nl; /* jump to next row of derivatives */ 3886 } 3887 PetscFunctionReturn(0); 3888 } 3889 3890 /* 3891 Special version for direct calls from Fortran 3892 */ 3893 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3894 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3895 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3896 #define matsetvaluesseqaij_ matsetvaluesseqaij 3897 #endif 3898 3899 /* Change these macros so can be used in void function */ 3900 #undef CHKERRQ 3901 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3902 #undef SETERRQ2 3903 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3904 3905 EXTERN_C_BEGIN 3906 #undef __FUNCT__ 3907 #define __FUNCT__ "matsetvaluesseqaij_" 3908 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3909 { 3910 Mat A = *AA; 3911 PetscInt m = *mm, n = *nn; 3912 InsertMode is = *isis; 3913 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3914 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3915 PetscInt *imax,*ai,*ailen; 3916 PetscErrorCode ierr; 3917 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3918 MatScalar *ap,value,*aa; 3919 PetscTruth ignorezeroentries = a->ignorezeroentries; 3920 PetscTruth roworiented = a->roworiented; 3921 3922 PetscFunctionBegin; 3923 ierr = MatPreallocated(A);CHKERRQ(ierr); 3924 imax = a->imax; 3925 ai = a->i; 3926 ailen = a->ilen; 3927 aj = a->j; 3928 aa = a->a; 3929 3930 for (k=0; k<m; k++) { /* loop over added rows */ 3931 row = im[k]; 3932 if (row < 0) continue; 3933 #if defined(PETSC_USE_DEBUG) 3934 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3935 #endif 3936 rp = aj + ai[row]; ap = aa + ai[row]; 3937 rmax = imax[row]; nrow = ailen[row]; 3938 low = 0; 3939 high = nrow; 3940 for (l=0; l<n; l++) { /* loop over added columns */ 3941 if (in[l] < 0) continue; 3942 #if defined(PETSC_USE_DEBUG) 3943 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3944 #endif 3945 col = in[l]; 3946 if (roworiented) { 3947 value = v[l + k*n]; 3948 } else { 3949 value = v[k + l*m]; 3950 } 3951 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3952 3953 if (col <= lastcol) low = 0; else high = nrow; 3954 lastcol = col; 3955 while (high-low > 5) { 3956 t = (low+high)/2; 3957 if (rp[t] > col) high = t; 3958 else low = t; 3959 } 3960 for (i=low; i<high; i++) { 3961 if (rp[i] > col) break; 3962 if (rp[i] == col) { 3963 if (is == ADD_VALUES) ap[i] += value; 3964 else ap[i] = value; 3965 goto noinsert; 3966 } 3967 } 3968 if (value == 0.0 && ignorezeroentries) goto noinsert; 3969 if (nonew == 1) goto noinsert; 3970 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3971 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3972 N = nrow++ - 1; a->nz++; high++; 3973 /* shift up all the later entries in this row */ 3974 for (ii=N; ii>=i; ii--) { 3975 rp[ii+1] = rp[ii]; 3976 ap[ii+1] = ap[ii]; 3977 } 3978 rp[i] = col; 3979 ap[i] = value; 3980 noinsert:; 3981 low = i + 1; 3982 } 3983 ailen[row] = nrow; 3984 } 3985 A->same_nonzero = PETSC_FALSE; 3986 PetscFunctionReturnVoid(); 3987 } 3988 EXTERN_C_END 3989