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