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