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 }; 2637 2638 EXTERN_C_BEGIN 2639 #undef __FUNCT__ 2640 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2641 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2642 { 2643 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2644 PetscInt i,nz,n; 2645 2646 PetscFunctionBegin; 2647 2648 nz = aij->maxnz; 2649 n = mat->rmap->n; 2650 for (i=0; i<nz; i++) { 2651 aij->j[i] = indices[i]; 2652 } 2653 aij->nz = nz; 2654 for (i=0; i<n; i++) { 2655 aij->ilen[i] = aij->imax[i]; 2656 } 2657 2658 PetscFunctionReturn(0); 2659 } 2660 EXTERN_C_END 2661 2662 #undef __FUNCT__ 2663 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2664 /*@ 2665 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2666 in the matrix. 2667 2668 Input Parameters: 2669 + mat - the SeqAIJ matrix 2670 - indices - the column indices 2671 2672 Level: advanced 2673 2674 Notes: 2675 This can be called if you have precomputed the nonzero structure of the 2676 matrix and want to provide it to the matrix object to improve the performance 2677 of the MatSetValues() operation. 2678 2679 You MUST have set the correct numbers of nonzeros per row in the call to 2680 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2681 2682 MUST be called before any calls to MatSetValues(); 2683 2684 The indices should start with zero, not one. 2685 2686 @*/ 2687 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2688 { 2689 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2690 2691 PetscFunctionBegin; 2692 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2693 PetscValidPointer(indices,2); 2694 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2695 if (f) { 2696 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2697 } else { 2698 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2699 } 2700 PetscFunctionReturn(0); 2701 } 2702 2703 /* ----------------------------------------------------------------------------------------*/ 2704 2705 EXTERN_C_BEGIN 2706 #undef __FUNCT__ 2707 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2708 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2709 { 2710 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2711 PetscErrorCode ierr; 2712 size_t nz = aij->i[mat->rmap->n]; 2713 2714 PetscFunctionBegin; 2715 if (aij->nonew != 1) { 2716 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2717 } 2718 2719 /* allocate space for values if not already there */ 2720 if (!aij->saved_values) { 2721 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2722 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2723 } 2724 2725 /* copy values over */ 2726 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2727 PetscFunctionReturn(0); 2728 } 2729 EXTERN_C_END 2730 2731 #undef __FUNCT__ 2732 #define __FUNCT__ "MatStoreValues" 2733 /*@ 2734 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2735 example, reuse of the linear part of a Jacobian, while recomputing the 2736 nonlinear portion. 2737 2738 Collect on Mat 2739 2740 Input Parameters: 2741 . mat - the matrix (currently only AIJ matrices support this option) 2742 2743 Level: advanced 2744 2745 Common Usage, with SNESSolve(): 2746 $ Create Jacobian matrix 2747 $ Set linear terms into matrix 2748 $ Apply boundary conditions to matrix, at this time matrix must have 2749 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2750 $ boundary conditions again will not change the nonzero structure 2751 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2752 $ ierr = MatStoreValues(mat); 2753 $ Call SNESSetJacobian() with matrix 2754 $ In your Jacobian routine 2755 $ ierr = MatRetrieveValues(mat); 2756 $ Set nonlinear terms in matrix 2757 2758 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2759 $ // build linear portion of Jacobian 2760 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2761 $ ierr = MatStoreValues(mat); 2762 $ loop over nonlinear iterations 2763 $ ierr = MatRetrieveValues(mat); 2764 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2765 $ // call MatAssemblyBegin/End() on matrix 2766 $ Solve linear system with Jacobian 2767 $ endloop 2768 2769 Notes: 2770 Matrix must already be assemblied before calling this routine 2771 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2772 calling this routine. 2773 2774 When this is called multiple times it overwrites the previous set of stored values 2775 and does not allocated additional space. 2776 2777 .seealso: MatRetrieveValues() 2778 2779 @*/ 2780 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2781 { 2782 PetscErrorCode ierr,(*f)(Mat); 2783 2784 PetscFunctionBegin; 2785 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2786 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2787 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2788 2789 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2790 if (f) { 2791 ierr = (*f)(mat);CHKERRQ(ierr); 2792 } else { 2793 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2794 } 2795 PetscFunctionReturn(0); 2796 } 2797 2798 EXTERN_C_BEGIN 2799 #undef __FUNCT__ 2800 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2801 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2802 { 2803 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2804 PetscErrorCode ierr; 2805 PetscInt nz = aij->i[mat->rmap->n]; 2806 2807 PetscFunctionBegin; 2808 if (aij->nonew != 1) { 2809 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2810 } 2811 if (!aij->saved_values) { 2812 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2813 } 2814 /* copy values over */ 2815 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2816 PetscFunctionReturn(0); 2817 } 2818 EXTERN_C_END 2819 2820 #undef __FUNCT__ 2821 #define __FUNCT__ "MatRetrieveValues" 2822 /*@ 2823 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2824 example, reuse of the linear part of a Jacobian, while recomputing the 2825 nonlinear portion. 2826 2827 Collect on Mat 2828 2829 Input Parameters: 2830 . mat - the matrix (currently on AIJ matrices support this option) 2831 2832 Level: advanced 2833 2834 .seealso: MatStoreValues() 2835 2836 @*/ 2837 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2838 { 2839 PetscErrorCode ierr,(*f)(Mat); 2840 2841 PetscFunctionBegin; 2842 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2843 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2844 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2845 2846 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2847 if (f) { 2848 ierr = (*f)(mat);CHKERRQ(ierr); 2849 } else { 2850 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2851 } 2852 PetscFunctionReturn(0); 2853 } 2854 2855 2856 /* --------------------------------------------------------------------------------*/ 2857 #undef __FUNCT__ 2858 #define __FUNCT__ "MatCreateSeqAIJ" 2859 /*@C 2860 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2861 (the default parallel PETSc format). For good matrix assembly performance 2862 the user should preallocate the matrix storage by setting the parameter nz 2863 (or the array nnz). By setting these parameters accurately, performance 2864 during matrix assembly can be increased by more than a factor of 50. 2865 2866 Collective on MPI_Comm 2867 2868 Input Parameters: 2869 + comm - MPI communicator, set to PETSC_COMM_SELF 2870 . m - number of rows 2871 . n - number of columns 2872 . nz - number of nonzeros per row (same for all rows) 2873 - nnz - array containing the number of nonzeros in the various rows 2874 (possibly different for each row) or PETSC_NULL 2875 2876 Output Parameter: 2877 . A - the matrix 2878 2879 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2880 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2881 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2882 2883 Notes: 2884 If nnz is given then nz is ignored 2885 2886 The AIJ format (also called the Yale sparse matrix format or 2887 compressed row storage), is fully compatible with standard Fortran 77 2888 storage. That is, the stored row and column indices can begin at 2889 either one (as in Fortran) or zero. See the users' manual for details. 2890 2891 Specify the preallocated storage with either nz or nnz (not both). 2892 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2893 allocation. For large problems you MUST preallocate memory or you 2894 will get TERRIBLE performance, see the users' manual chapter on matrices. 2895 2896 By default, this format uses inodes (identical nodes) when possible, to 2897 improve numerical efficiency of matrix-vector products and solves. We 2898 search for consecutive rows with the same nonzero structure, thereby 2899 reusing matrix information to achieve increased efficiency. 2900 2901 Options Database Keys: 2902 + -mat_no_inode - Do not use inodes 2903 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2904 - -mat_aij_oneindex - Internally use indexing starting at 1 2905 rather than 0. Note that when calling MatSetValues(), 2906 the user still MUST index entries starting at 0! 2907 2908 Level: intermediate 2909 2910 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2911 2912 @*/ 2913 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2914 { 2915 PetscErrorCode ierr; 2916 2917 PetscFunctionBegin; 2918 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2919 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2920 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2921 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2922 PetscFunctionReturn(0); 2923 } 2924 2925 #undef __FUNCT__ 2926 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2927 /*@C 2928 MatSeqAIJSetPreallocation - For good matrix assembly performance 2929 the user should preallocate the matrix storage by setting the parameter nz 2930 (or the array nnz). By setting these parameters accurately, performance 2931 during matrix assembly can be increased by more than a factor of 50. 2932 2933 Collective on MPI_Comm 2934 2935 Input Parameters: 2936 + B - The matrix-free 2937 . nz - number of nonzeros per row (same for all rows) 2938 - nnz - array containing the number of nonzeros in the various rows 2939 (possibly different for each row) or PETSC_NULL 2940 2941 Notes: 2942 If nnz is given then nz is ignored 2943 2944 The AIJ format (also called the Yale sparse matrix format or 2945 compressed row storage), is fully compatible with standard Fortran 77 2946 storage. That is, the stored row and column indices can begin at 2947 either one (as in Fortran) or zero. See the users' manual for details. 2948 2949 Specify the preallocated storage with either nz or nnz (not both). 2950 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2951 allocation. For large problems you MUST preallocate memory or you 2952 will get TERRIBLE performance, see the users' manual chapter on matrices. 2953 2954 You can call MatGetInfo() to get information on how effective the preallocation was; 2955 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2956 You can also run with the option -info and look for messages with the string 2957 malloc in them to see if additional memory allocation was needed. 2958 2959 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2960 entries or columns indices 2961 2962 By default, this format uses inodes (identical nodes) when possible, to 2963 improve numerical efficiency of matrix-vector products and solves. We 2964 search for consecutive rows with the same nonzero structure, thereby 2965 reusing matrix information to achieve increased efficiency. 2966 2967 Options Database Keys: 2968 + -mat_no_inode - Do not use inodes 2969 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2970 - -mat_aij_oneindex - Internally use indexing starting at 1 2971 rather than 0. Note that when calling MatSetValues(), 2972 the user still MUST index entries starting at 0! 2973 2974 Level: intermediate 2975 2976 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2977 2978 @*/ 2979 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2980 { 2981 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2982 2983 PetscFunctionBegin; 2984 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2985 if (f) { 2986 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2987 } 2988 PetscFunctionReturn(0); 2989 } 2990 2991 EXTERN_C_BEGIN 2992 #undef __FUNCT__ 2993 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2994 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2995 { 2996 Mat_SeqAIJ *b; 2997 PetscTruth skipallocation = PETSC_FALSE; 2998 PetscErrorCode ierr; 2999 PetscInt i; 3000 3001 PetscFunctionBegin; 3002 3003 if (nz == MAT_SKIP_ALLOCATION) { 3004 skipallocation = PETSC_TRUE; 3005 nz = 0; 3006 } 3007 3008 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 3009 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 3010 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3011 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3012 3013 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3014 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3015 if (nnz) { 3016 for (i=0; i<B->rmap->n; i++) { 3017 if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 3018 if (nnz[i] > B->cmap->n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->cmap->n); 3019 } 3020 } 3021 3022 B->preallocated = PETSC_TRUE; 3023 b = (Mat_SeqAIJ*)B->data; 3024 3025 if (!skipallocation) { 3026 if (!b->imax) { 3027 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3028 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3029 } 3030 if (!nnz) { 3031 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3032 else if (nz <= 0) nz = 1; 3033 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3034 nz = nz*B->rmap->n; 3035 } else { 3036 nz = 0; 3037 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3038 } 3039 /* b->ilen will count nonzeros in each row so far. */ 3040 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3041 3042 /* allocate the matrix space */ 3043 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3044 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3045 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3046 b->i[0] = 0; 3047 for (i=1; i<B->rmap->n+1; i++) { 3048 b->i[i] = b->i[i-1] + b->imax[i-1]; 3049 } 3050 b->singlemalloc = PETSC_TRUE; 3051 b->free_a = PETSC_TRUE; 3052 b->free_ij = PETSC_TRUE; 3053 } else { 3054 b->free_a = PETSC_FALSE; 3055 b->free_ij = PETSC_FALSE; 3056 } 3057 3058 b->nz = 0; 3059 b->maxnz = nz; 3060 B->info.nz_unneeded = (double)b->maxnz; 3061 PetscFunctionReturn(0); 3062 } 3063 EXTERN_C_END 3064 3065 #undef __FUNCT__ 3066 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3067 /*@ 3068 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3069 3070 Input Parameters: 3071 + B - the matrix 3072 . i - the indices into j for the start of each row (starts with zero) 3073 . j - the column indices for each row (starts with zero) these must be sorted for each row 3074 - v - optional values in the matrix 3075 3076 Level: developer 3077 3078 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3079 3080 .keywords: matrix, aij, compressed row, sparse, sequential 3081 3082 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3083 @*/ 3084 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3085 { 3086 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3087 PetscErrorCode ierr; 3088 3089 PetscFunctionBegin; 3090 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3091 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3092 if (f) { 3093 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3094 } 3095 PetscFunctionReturn(0); 3096 } 3097 3098 EXTERN_C_BEGIN 3099 #undef __FUNCT__ 3100 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3101 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3102 { 3103 PetscInt i; 3104 PetscInt m,n; 3105 PetscInt nz; 3106 PetscInt *nnz, nz_max = 0; 3107 PetscScalar *values; 3108 PetscErrorCode ierr; 3109 3110 PetscFunctionBegin; 3111 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3112 3113 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3114 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3115 for(i = 0; i < m; i++) { 3116 nz = Ii[i+1]- Ii[i]; 3117 nz_max = PetscMax(nz_max, nz); 3118 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3119 nnz[i] = nz; 3120 } 3121 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3122 ierr = PetscFree(nnz);CHKERRQ(ierr); 3123 3124 if (v) { 3125 values = (PetscScalar*) v; 3126 } else { 3127 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3128 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3129 } 3130 3131 for(i = 0; i < m; i++) { 3132 nz = Ii[i+1] - Ii[i]; 3133 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3134 } 3135 3136 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3137 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3138 3139 if (!v) { 3140 ierr = PetscFree(values);CHKERRQ(ierr); 3141 } 3142 PetscFunctionReturn(0); 3143 } 3144 EXTERN_C_END 3145 3146 #include "../src/mat/impls/dense/seq/dense.h" 3147 #include "private/petscaxpy.h" 3148 3149 #undef __FUNCT__ 3150 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3151 /* 3152 Computes (B'*A')' since computing B*A directly is untenable 3153 3154 n p p 3155 ( ) ( ) ( ) 3156 m ( A ) * n ( B ) = m ( C ) 3157 ( ) ( ) ( ) 3158 3159 */ 3160 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3161 { 3162 PetscErrorCode ierr; 3163 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3164 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3165 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3166 PetscInt i,n,m,q,p; 3167 const PetscInt *ii,*idx; 3168 const PetscScalar *b,*a,*a_q; 3169 PetscScalar *c,*c_q; 3170 3171 PetscFunctionBegin; 3172 m = A->rmap->n; 3173 n = A->cmap->n; 3174 p = B->cmap->n; 3175 a = sub_a->v; 3176 b = sub_b->a; 3177 c = sub_c->v; 3178 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3179 3180 ii = sub_b->i; 3181 idx = sub_b->j; 3182 for (i=0; i<n; i++) { 3183 q = ii[i+1] - ii[i]; 3184 while (q-->0) { 3185 c_q = c + m*(*idx); 3186 a_q = a + m*i; 3187 PetscAXPY(c_q,*b,a_q,m); 3188 idx++; 3189 b++; 3190 } 3191 } 3192 PetscFunctionReturn(0); 3193 } 3194 3195 #undef __FUNCT__ 3196 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3197 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3198 { 3199 PetscErrorCode ierr; 3200 PetscInt m=A->rmap->n,n=B->cmap->n; 3201 Mat Cmat; 3202 3203 PetscFunctionBegin; 3204 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 3205 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3206 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3207 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3208 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3209 Cmat->assembled = PETSC_TRUE; 3210 *C = Cmat; 3211 PetscFunctionReturn(0); 3212 } 3213 3214 /* ----------------------------------------------------------------*/ 3215 #undef __FUNCT__ 3216 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3217 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3218 { 3219 PetscErrorCode ierr; 3220 3221 PetscFunctionBegin; 3222 if (scall == MAT_INITIAL_MATRIX){ 3223 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3224 } 3225 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3226 PetscFunctionReturn(0); 3227 } 3228 3229 3230 /*MC 3231 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3232 based on compressed sparse row format. 3233 3234 Options Database Keys: 3235 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3236 3237 Level: beginner 3238 3239 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3240 M*/ 3241 3242 EXTERN_C_BEGIN 3243 #if defined(PETSC_HAVE_PASTIX) 3244 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3245 #endif 3246 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3247 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3248 #endif 3249 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3250 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3251 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3252 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3253 #if defined(PETSC_HAVE_MUMPS) 3254 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3255 #endif 3256 #if defined(PETSC_HAVE_SUPERLU) 3257 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3258 #endif 3259 #if defined(PETSC_HAVE_SUPERLU_DIST) 3260 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3261 #endif 3262 #if defined(PETSC_HAVE_SPOOLES) 3263 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3264 #endif 3265 #if defined(PETSC_HAVE_UMFPACK) 3266 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3267 #endif 3268 #if defined(PETSC_HAVE_CHOLMOD) 3269 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3270 #endif 3271 #if defined(PETSC_HAVE_LUSOL) 3272 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3273 #endif 3274 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3275 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3276 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEnginePut_SeqAIJ(PetscObject,void*); 3277 extern PetscErrorCode PETSCMAT_DLLEXPORT MatlabEngineGet_SeqAIJ(PetscObject,void*); 3278 #endif 3279 EXTERN_C_END 3280 3281 3282 EXTERN_C_BEGIN 3283 #undef __FUNCT__ 3284 #define __FUNCT__ "MatCreate_SeqAIJ" 3285 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3286 { 3287 Mat_SeqAIJ *b; 3288 PetscErrorCode ierr; 3289 PetscMPIInt size; 3290 3291 PetscFunctionBegin; 3292 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3293 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3294 3295 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3296 B->data = (void*)b; 3297 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3298 B->mapping = 0; 3299 b->row = 0; 3300 b->col = 0; 3301 b->icol = 0; 3302 b->reallocs = 0; 3303 b->ignorezeroentries = PETSC_FALSE; 3304 b->roworiented = PETSC_TRUE; 3305 b->nonew = 0; 3306 b->diag = 0; 3307 b->solve_work = 0; 3308 B->spptr = 0; 3309 b->saved_values = 0; 3310 b->idiag = 0; 3311 b->mdiag = 0; 3312 b->ssor_work = 0; 3313 b->omega = 1.0; 3314 b->fshift = 0.0; 3315 b->idiagvalid = PETSC_FALSE; 3316 b->keepnonzeropattern = PETSC_FALSE; 3317 b->xtoy = 0; 3318 b->XtoY = 0; 3319 b->compressedrow.use = PETSC_FALSE; 3320 b->compressedrow.nrows = B->rmap->n; 3321 b->compressedrow.i = PETSC_NULL; 3322 b->compressedrow.rindex = PETSC_NULL; 3323 b->compressedrow.checked = PETSC_FALSE; 3324 B->same_nonzero = PETSC_FALSE; 3325 3326 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3327 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3328 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C", 3329 "MatGetFactor_seqaij_matlab", 3330 MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3331 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3332 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3333 #endif 3334 #if defined(PETSC_HAVE_PASTIX) 3335 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3336 "MatGetFactor_seqaij_pastix", 3337 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3338 #endif 3339 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3340 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3341 "MatGetFactor_seqaij_essl", 3342 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3343 #endif 3344 #if defined(PETSC_HAVE_SUPERLU) 3345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3346 "MatGetFactor_seqaij_superlu", 3347 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3348 #endif 3349 #if defined(PETSC_HAVE_SUPERLU_DIST) 3350 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3351 "MatGetFactor_seqaij_superlu_dist", 3352 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3353 #endif 3354 #if defined(PETSC_HAVE_SPOOLES) 3355 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3356 "MatGetFactor_seqaij_spooles", 3357 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3358 #endif 3359 #if defined(PETSC_HAVE_MUMPS) 3360 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3361 "MatGetFactor_aij_mumps", 3362 MatGetFactor_aij_mumps);CHKERRQ(ierr); 3363 #endif 3364 #if defined(PETSC_HAVE_UMFPACK) 3365 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3366 "MatGetFactor_seqaij_umfpack", 3367 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3368 #endif 3369 #if defined(PETSC_HAVE_CHOLMOD) 3370 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 3371 "MatGetFactor_seqaij_cholmod", 3372 MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3373 #endif 3374 #if defined(PETSC_HAVE_LUSOL) 3375 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3376 "MatGetFactor_seqaij_lusol", 3377 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3378 #endif 3379 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3380 "MatGetFactor_seqaij_petsc", 3381 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3382 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3383 "MatGetFactorAvailable_seqaij_petsc", 3384 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3385 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3386 "MatGetFactor_seqaij_bas", 3387 MatGetFactor_seqaij_bas); 3388 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3389 "MatSeqAIJSetColumnIndices_SeqAIJ", 3390 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3391 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3392 "MatStoreValues_SeqAIJ", 3393 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3394 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3395 "MatRetrieveValues_SeqAIJ", 3396 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3397 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3398 "MatConvert_SeqAIJ_SeqSBAIJ", 3399 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3400 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3401 "MatConvert_SeqAIJ_SeqBAIJ", 3402 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3403 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3404 "MatConvert_SeqAIJ_SeqCSRPERM", 3405 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3406 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3407 "MatConvert_SeqAIJ_SeqCRL", 3408 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3409 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3410 "MatIsTranspose_SeqAIJ", 3411 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3412 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3413 "MatIsHermitianTranspose_SeqAIJ", 3414 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3415 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3416 "MatSeqAIJSetPreallocation_SeqAIJ", 3417 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3418 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3419 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3420 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3421 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3422 "MatReorderForNonzeroDiagonal_SeqAIJ", 3423 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3424 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3425 "MatMatMult_SeqDense_SeqAIJ", 3426 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3427 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3428 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3429 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3430 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3431 "MatMatMultNumeric_SeqDense_SeqAIJ", 3432 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3433 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3434 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3435 PetscFunctionReturn(0); 3436 } 3437 EXTERN_C_END 3438 3439 #undef __FUNCT__ 3440 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3441 /* 3442 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3443 */ 3444 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3445 { 3446 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3447 PetscErrorCode ierr; 3448 PetscInt i,m = A->rmap->n; 3449 3450 PetscFunctionBegin; 3451 c = (Mat_SeqAIJ*)C->data; 3452 3453 C->factortype = A->factortype; 3454 c->row = 0; 3455 c->col = 0; 3456 c->icol = 0; 3457 c->reallocs = 0; 3458 3459 C->assembled = PETSC_TRUE; 3460 3461 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3462 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3463 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3464 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3465 3466 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3467 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3468 for (i=0; i<m; i++) { 3469 c->imax[i] = a->imax[i]; 3470 c->ilen[i] = a->ilen[i]; 3471 } 3472 3473 /* allocate the matrix space */ 3474 if (mallocmatspace){ 3475 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3476 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3477 c->singlemalloc = PETSC_TRUE; 3478 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3479 if (m > 0) { 3480 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3481 if (cpvalues == MAT_COPY_VALUES) { 3482 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3483 } else { 3484 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3485 } 3486 } 3487 } 3488 3489 c->ignorezeroentries = a->ignorezeroentries; 3490 c->roworiented = a->roworiented; 3491 c->nonew = a->nonew; 3492 if (a->diag) { 3493 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3494 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3495 for (i=0; i<m; i++) { 3496 c->diag[i] = a->diag[i]; 3497 } 3498 } else c->diag = 0; 3499 c->solve_work = 0; 3500 c->saved_values = 0; 3501 c->idiag = 0; 3502 c->ssor_work = 0; 3503 c->keepnonzeropattern = a->keepnonzeropattern; 3504 c->free_a = PETSC_TRUE; 3505 c->free_ij = PETSC_TRUE; 3506 c->xtoy = 0; 3507 c->XtoY = 0; 3508 3509 c->nz = a->nz; 3510 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 3511 C->preallocated = PETSC_TRUE; 3512 3513 c->compressedrow.use = a->compressedrow.use; 3514 c->compressedrow.nrows = a->compressedrow.nrows; 3515 c->compressedrow.checked = a->compressedrow.checked; 3516 if (a->compressedrow.checked && a->compressedrow.use){ 3517 i = a->compressedrow.nrows; 3518 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3519 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3520 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3521 } else { 3522 c->compressedrow.use = PETSC_FALSE; 3523 c->compressedrow.i = PETSC_NULL; 3524 c->compressedrow.rindex = PETSC_NULL; 3525 } 3526 C->same_nonzero = A->same_nonzero; 3527 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3528 3529 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3530 PetscFunctionReturn(0); 3531 } 3532 3533 #undef __FUNCT__ 3534 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3535 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3536 { 3537 PetscErrorCode ierr; 3538 3539 PetscFunctionBegin; 3540 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3541 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3542 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3543 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3544 PetscFunctionReturn(0); 3545 } 3546 3547 #undef __FUNCT__ 3548 #define __FUNCT__ "MatLoad_SeqAIJ" 3549 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer,Mat newMat) 3550 { 3551 Mat_SeqAIJ *a; 3552 PetscErrorCode ierr; 3553 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 3554 int fd; 3555 PetscMPIInt size; 3556 MPI_Comm comm; 3557 3558 PetscFunctionBegin; 3559 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3560 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3561 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 3562 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3563 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3564 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3565 M = header[1]; N = header[2]; nz = header[3]; 3566 3567 if (nz < 0) { 3568 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3569 } 3570 3571 /* read in row lengths */ 3572 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3573 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3574 3575 /* check if sum of rowlengths is same as nz */ 3576 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3577 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); 3578 3579 /* set global size if not set already*/ 3580 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 3581 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3582 } else { 3583 /* if sizes and type are already set, check if the vector global sizes are correct */ 3584 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 3585 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); 3586 } 3587 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 3588 a = (Mat_SeqAIJ*)newMat->data; 3589 3590 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3591 3592 /* read in nonzero values */ 3593 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3594 3595 /* set matrix "i" values */ 3596 a->i[0] = 0; 3597 for (i=1; i<= M; i++) { 3598 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3599 a->ilen[i-1] = rowlengths[i-1]; 3600 } 3601 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3602 3603 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3604 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3605 PetscFunctionReturn(0); 3606 } 3607 3608 #undef __FUNCT__ 3609 #define __FUNCT__ "MatEqual_SeqAIJ" 3610 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3611 { 3612 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3613 PetscErrorCode ierr; 3614 #if defined(PETSC_USE_COMPLEX) 3615 PetscInt k; 3616 #endif 3617 3618 PetscFunctionBegin; 3619 /* If the matrix dimensions are not equal,or no of nonzeros */ 3620 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3621 *flg = PETSC_FALSE; 3622 PetscFunctionReturn(0); 3623 } 3624 3625 /* if the a->i are the same */ 3626 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3627 if (!*flg) PetscFunctionReturn(0); 3628 3629 /* if a->j are the same */ 3630 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3631 if (!*flg) PetscFunctionReturn(0); 3632 3633 /* if a->a are the same */ 3634 #if defined(PETSC_USE_COMPLEX) 3635 for (k=0; k<a->nz; k++){ 3636 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3637 *flg = PETSC_FALSE; 3638 PetscFunctionReturn(0); 3639 } 3640 } 3641 #else 3642 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3643 #endif 3644 PetscFunctionReturn(0); 3645 } 3646 3647 #undef __FUNCT__ 3648 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3649 /*@ 3650 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3651 provided by the user. 3652 3653 Collective on MPI_Comm 3654 3655 Input Parameters: 3656 + comm - must be an MPI communicator of size 1 3657 . m - number of rows 3658 . n - number of columns 3659 . i - row indices 3660 . j - column indices 3661 - a - matrix values 3662 3663 Output Parameter: 3664 . mat - the matrix 3665 3666 Level: intermediate 3667 3668 Notes: 3669 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3670 once the matrix is destroyed 3671 3672 You cannot set new nonzero locations into this matrix, that will generate an error. 3673 3674 The i and j indices are 0 based 3675 3676 The format which is used for the sparse matrix input, is equivalent to a 3677 row-major ordering.. i.e for the following matrix, the input data expected is 3678 as shown: 3679 3680 1 0 0 3681 2 0 3 3682 4 5 6 3683 3684 i = {0,1,3,6} [size = nrow+1 = 3+1] 3685 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3686 v = {1,2,3,4,5,6} [size = nz = 6] 3687 3688 3689 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3690 3691 @*/ 3692 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3693 { 3694 PetscErrorCode ierr; 3695 PetscInt ii; 3696 Mat_SeqAIJ *aij; 3697 #if defined(PETSC_USE_DEBUG) 3698 PetscInt jj; 3699 #endif 3700 3701 PetscFunctionBegin; 3702 if (i[0]) { 3703 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3704 } 3705 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3706 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3707 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3708 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3709 aij = (Mat_SeqAIJ*)(*mat)->data; 3710 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3711 3712 aij->i = i; 3713 aij->j = j; 3714 aij->a = a; 3715 aij->singlemalloc = PETSC_FALSE; 3716 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3717 aij->free_a = PETSC_FALSE; 3718 aij->free_ij = PETSC_FALSE; 3719 3720 for (ii=0; ii<m; ii++) { 3721 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3722 #if defined(PETSC_USE_DEBUG) 3723 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]); 3724 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3725 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); 3726 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); 3727 } 3728 #endif 3729 } 3730 #if defined(PETSC_USE_DEBUG) 3731 for (ii=0; ii<aij->i[m]; ii++) { 3732 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3733 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]); 3734 } 3735 #endif 3736 3737 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3738 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3739 PetscFunctionReturn(0); 3740 } 3741 3742 #undef __FUNCT__ 3743 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3744 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3745 { 3746 PetscErrorCode ierr; 3747 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3748 3749 PetscFunctionBegin; 3750 if (coloring->ctype == IS_COLORING_GLOBAL) { 3751 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3752 a->coloring = coloring; 3753 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3754 PetscInt i,*larray; 3755 ISColoring ocoloring; 3756 ISColoringValue *colors; 3757 3758 /* set coloring for diagonal portion */ 3759 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3760 for (i=0; i<A->cmap->n; i++) { 3761 larray[i] = i; 3762 } 3763 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3764 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3765 for (i=0; i<A->cmap->n; i++) { 3766 colors[i] = coloring->colors[larray[i]]; 3767 } 3768 ierr = PetscFree(larray);CHKERRQ(ierr); 3769 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3770 a->coloring = ocoloring; 3771 } 3772 PetscFunctionReturn(0); 3773 } 3774 3775 #if defined(PETSC_HAVE_ADIC) 3776 EXTERN_C_BEGIN 3777 #include "adic/ad_utils.h" 3778 EXTERN_C_END 3779 3780 #undef __FUNCT__ 3781 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3782 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3783 { 3784 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3785 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3786 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3787 ISColoringValue *color; 3788 3789 PetscFunctionBegin; 3790 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3791 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3792 color = a->coloring->colors; 3793 /* loop over rows */ 3794 for (i=0; i<m; i++) { 3795 nz = ii[i+1] - ii[i]; 3796 /* loop over columns putting computed value into matrix */ 3797 for (j=0; j<nz; j++) { 3798 *v++ = values[color[*jj++]]; 3799 } 3800 values += nlen; /* jump to next row of derivatives */ 3801 } 3802 PetscFunctionReturn(0); 3803 } 3804 #endif 3805 3806 #undef __FUNCT__ 3807 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3808 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3809 { 3810 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3811 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3812 MatScalar *v = a->a; 3813 PetscScalar *values = (PetscScalar *)advalues; 3814 ISColoringValue *color; 3815 3816 PetscFunctionBegin; 3817 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3818 color = a->coloring->colors; 3819 /* loop over rows */ 3820 for (i=0; i<m; i++) { 3821 nz = ii[i+1] - ii[i]; 3822 /* loop over columns putting computed value into matrix */ 3823 for (j=0; j<nz; j++) { 3824 *v++ = values[color[*jj++]]; 3825 } 3826 values += nl; /* jump to next row of derivatives */ 3827 } 3828 PetscFunctionReturn(0); 3829 } 3830 3831 /* 3832 Special version for direct calls from Fortran 3833 */ 3834 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3835 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3836 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3837 #define matsetvaluesseqaij_ matsetvaluesseqaij 3838 #endif 3839 3840 /* Change these macros so can be used in void function */ 3841 #undef CHKERRQ 3842 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3843 #undef SETERRQ2 3844 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 3845 3846 EXTERN_C_BEGIN 3847 #undef __FUNCT__ 3848 #define __FUNCT__ "matsetvaluesseqaij_" 3849 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3850 { 3851 Mat A = *AA; 3852 PetscInt m = *mm, n = *nn; 3853 InsertMode is = *isis; 3854 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3855 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3856 PetscInt *imax,*ai,*ailen; 3857 PetscErrorCode ierr; 3858 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3859 MatScalar *ap,value,*aa; 3860 PetscTruth ignorezeroentries = a->ignorezeroentries; 3861 PetscTruth roworiented = a->roworiented; 3862 3863 PetscFunctionBegin; 3864 ierr = MatPreallocated(A);CHKERRQ(ierr); 3865 imax = a->imax; 3866 ai = a->i; 3867 ailen = a->ilen; 3868 aj = a->j; 3869 aa = a->a; 3870 3871 for (k=0; k<m; k++) { /* loop over added rows */ 3872 row = im[k]; 3873 if (row < 0) continue; 3874 #if defined(PETSC_USE_DEBUG) 3875 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3876 #endif 3877 rp = aj + ai[row]; ap = aa + ai[row]; 3878 rmax = imax[row]; nrow = ailen[row]; 3879 low = 0; 3880 high = nrow; 3881 for (l=0; l<n; l++) { /* loop over added columns */ 3882 if (in[l] < 0) continue; 3883 #if defined(PETSC_USE_DEBUG) 3884 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3885 #endif 3886 col = in[l]; 3887 if (roworiented) { 3888 value = v[l + k*n]; 3889 } else { 3890 value = v[k + l*m]; 3891 } 3892 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3893 3894 if (col <= lastcol) low = 0; else high = nrow; 3895 lastcol = col; 3896 while (high-low > 5) { 3897 t = (low+high)/2; 3898 if (rp[t] > col) high = t; 3899 else low = t; 3900 } 3901 for (i=low; i<high; i++) { 3902 if (rp[i] > col) break; 3903 if (rp[i] == col) { 3904 if (is == ADD_VALUES) ap[i] += value; 3905 else ap[i] = value; 3906 goto noinsert; 3907 } 3908 } 3909 if (value == 0.0 && ignorezeroentries) goto noinsert; 3910 if (nonew == 1) goto noinsert; 3911 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3912 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3913 N = nrow++ - 1; a->nz++; high++; 3914 /* shift up all the later entries in this row */ 3915 for (ii=N; ii>=i; ii--) { 3916 rp[ii+1] = rp[ii]; 3917 ap[ii+1] = ap[ii]; 3918 } 3919 rp[i] = col; 3920 ap[i] = value; 3921 noinsert:; 3922 low = i + 1; 3923 } 3924 ailen[row] = nrow; 3925 } 3926 A->same_nonzero = PETSC_FALSE; 3927 PetscFunctionReturnVoid(); 3928 } 3929 EXTERN_C_END 3930