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