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