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