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