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 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 1073 1074 diag = a->diag; 1075 if (!a->idiag) { 1076 ierr = PetscMalloc(3*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1077 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1078 a->ssor = a->idiag + m; 1079 mdiag = a->ssor + m; 1080 v = a->a; 1081 1082 /* this is wrong when fshift omega changes each iteration */ 1083 if (omega == 1.0 && !fshift) { 1084 for (i=0; i<m; i++) { 1085 mdiag[i] = v[diag[i]]; 1086 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1087 a->idiag[i] = 1.0/v[diag[i]]; 1088 } 1089 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1090 } else { 1091 for (i=0; i<m; i++) { 1092 mdiag[i] = v[diag[i]]; 1093 a->idiag[i] = omega/(fshift + v[diag[i]]); 1094 } 1095 ierr = PetscLogFlops(2*m);CHKERRQ(ierr); 1096 } 1097 } 1098 t = a->ssor; 1099 idiag = a->idiag; 1100 mdiag = a->idiag + 2*m; 1101 1102 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1103 if (xx != bb) { 1104 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1105 } else { 1106 b = x; 1107 } 1108 1109 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1110 xs = x; 1111 if (flag == SOR_APPLY_UPPER) { 1112 /* apply (U + D/omega) to the vector */ 1113 bs = b; 1114 for (i=0; i<m; i++) { 1115 d = fshift + a->a[diag[i]]; 1116 n = a->i[i+1] - diag[i] - 1; 1117 idx = a->j + diag[i] + 1; 1118 v = a->a + diag[i] + 1; 1119 sum = b[i]*d/omega; 1120 SPARSEDENSEDOT(sum,bs,v,idx,n); 1121 x[i] = sum; 1122 } 1123 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1124 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1125 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1126 PetscFunctionReturn(0); 1127 } 1128 1129 1130 /* Let A = L + U + D; where L is lower trianglar, 1131 U is upper triangular, E is diagonal; This routine applies 1132 1133 (L + E)^{-1} A (U + E)^{-1} 1134 1135 to a vector efficiently using Eisenstat's trick. This is for 1136 the case of SSOR preconditioner, so E is D/omega where omega 1137 is the relaxation factor. 1138 */ 1139 1140 if (flag == SOR_APPLY_LOWER) { 1141 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1142 } else if (flag & SOR_EISENSTAT) { 1143 /* Let A = L + U + D; where L is lower trianglar, 1144 U is upper triangular, E is diagonal; This routine applies 1145 1146 (L + E)^{-1} A (U + E)^{-1} 1147 1148 to a vector efficiently using Eisenstat's trick. This is for 1149 the case of SSOR preconditioner, so E is D/omega where omega 1150 is the relaxation factor. 1151 */ 1152 scale = (2.0/omega) - 1.0; 1153 1154 /* x = (E + U)^{-1} b */ 1155 for (i=m-1; i>=0; i--) { 1156 n = a->i[i+1] - diag[i] - 1; 1157 idx = a->j + diag[i] + 1; 1158 v = a->a + diag[i] + 1; 1159 sum = b[i]; 1160 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1161 x[i] = sum*idiag[i]; 1162 } 1163 1164 /* t = b - (2*E - D)x */ 1165 v = a->a; 1166 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1167 1168 /* t = (E + L)^{-1}t */ 1169 ts = t; 1170 diag = a->diag; 1171 for (i=0; i<m; i++) { 1172 n = diag[i] - a->i[i]; 1173 idx = a->j + a->i[i]; 1174 v = a->a + a->i[i]; 1175 sum = t[i]; 1176 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1177 t[i] = sum*idiag[i]; 1178 /* x = x + t */ 1179 x[i] += t[i]; 1180 } 1181 1182 ierr = PetscLogFlops(6*m-1 + 2*a->nz);CHKERRQ(ierr); 1183 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1184 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1185 PetscFunctionReturn(0); 1186 } 1187 if (flag & SOR_ZERO_INITIAL_GUESS) { 1188 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1189 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1190 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)b); 1191 #else 1192 for (i=0; i<m; i++) { 1193 n = diag[i] - a->i[i]; 1194 idx = a->j + a->i[i]; 1195 v = a->a + a->i[i]; 1196 sum = b[i]; 1197 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1198 x[i] = sum*idiag[i]; 1199 } 1200 #endif 1201 xb = x; 1202 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1203 } else xb = b; 1204 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1205 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1206 for (i=0; i<m; i++) { 1207 x[i] *= mdiag[i]; 1208 } 1209 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1210 } 1211 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1212 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1213 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,idiag,a->a,(void*)xb); 1214 #else 1215 for (i=m-1; i>=0; i--) { 1216 n = a->i[i+1] - diag[i] - 1; 1217 idx = a->j + diag[i] + 1; 1218 v = a->a + diag[i] + 1; 1219 sum = xb[i]; 1220 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1221 x[i] = sum*idiag[i]; 1222 } 1223 #endif 1224 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1225 } 1226 its--; 1227 } 1228 while (its--) { 1229 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1230 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1231 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1232 #else 1233 for (i=0; i<m; i++) { 1234 n = a->i[i+1] - a->i[i]; 1235 idx = a->j + a->i[i]; 1236 v = a->a + a->i[i]; 1237 sum = b[i]; 1238 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1239 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1240 } 1241 #endif 1242 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1243 } 1244 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1245 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1246 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,(PetscInt*)diag,a->a,(void*)b); 1247 #else 1248 for (i=m-1; i>=0; i--) { 1249 n = a->i[i+1] - a->i[i]; 1250 idx = a->j + a->i[i]; 1251 v = a->a + a->i[i]; 1252 sum = b[i]; 1253 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1254 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1255 } 1256 #endif 1257 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1258 } 1259 } 1260 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1261 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1262 PetscFunctionReturn(0); 1263 } 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1267 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1268 { 1269 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1270 1271 PetscFunctionBegin; 1272 info->rows_global = (double)A->rmap.n; 1273 info->columns_global = (double)A->cmap.n; 1274 info->rows_local = (double)A->rmap.n; 1275 info->columns_local = (double)A->cmap.n; 1276 info->block_size = 1.0; 1277 info->nz_allocated = (double)a->maxnz; 1278 info->nz_used = (double)a->nz; 1279 info->nz_unneeded = (double)(a->maxnz - a->nz); 1280 info->assemblies = (double)A->num_ass; 1281 info->mallocs = (double)a->reallocs; 1282 info->memory = A->mem; 1283 if (A->factor) { 1284 info->fill_ratio_given = A->info.fill_ratio_given; 1285 info->fill_ratio_needed = A->info.fill_ratio_needed; 1286 info->factor_mallocs = A->info.factor_mallocs; 1287 } else { 1288 info->fill_ratio_given = 0; 1289 info->fill_ratio_needed = 0; 1290 info->factor_mallocs = 0; 1291 } 1292 PetscFunctionReturn(0); 1293 } 1294 1295 #undef __FUNCT__ 1296 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1297 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1298 { 1299 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1300 PetscInt i,m = A->rmap.n - 1,d; 1301 PetscErrorCode ierr; 1302 PetscTruth missing; 1303 1304 PetscFunctionBegin; 1305 if (a->keepzeroedrows) { 1306 for (i=0; i<N; i++) { 1307 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1308 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1309 } 1310 if (diag != 0.0) { 1311 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1312 if (missing) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1313 for (i=0; i<N; i++) { 1314 a->a[a->diag[rows[i]]] = diag; 1315 } 1316 } 1317 A->same_nonzero = PETSC_TRUE; 1318 } else { 1319 if (diag != 0.0) { 1320 for (i=0; i<N; i++) { 1321 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1322 if (a->ilen[rows[i]] > 0) { 1323 a->ilen[rows[i]] = 1; 1324 a->a[a->i[rows[i]]] = diag; 1325 a->j[a->i[rows[i]]] = rows[i]; 1326 } else { /* in case row was completely empty */ 1327 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1328 } 1329 } 1330 } else { 1331 for (i=0; i<N; i++) { 1332 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1333 a->ilen[rows[i]] = 0; 1334 } 1335 } 1336 A->same_nonzero = PETSC_FALSE; 1337 } 1338 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1339 PetscFunctionReturn(0); 1340 } 1341 1342 #undef __FUNCT__ 1343 #define __FUNCT__ "MatGetRow_SeqAIJ" 1344 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1345 { 1346 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1347 PetscInt *itmp; 1348 1349 PetscFunctionBegin; 1350 if (row < 0 || row >= A->rmap.n) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1351 1352 *nz = a->i[row+1] - a->i[row]; 1353 if (v) *v = a->a + a->i[row]; 1354 if (idx) { 1355 itmp = a->j + a->i[row]; 1356 if (*nz) { 1357 *idx = itmp; 1358 } 1359 else *idx = 0; 1360 } 1361 PetscFunctionReturn(0); 1362 } 1363 1364 /* remove this function? */ 1365 #undef __FUNCT__ 1366 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1367 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1368 { 1369 PetscFunctionBegin; 1370 PetscFunctionReturn(0); 1371 } 1372 1373 #undef __FUNCT__ 1374 #define __FUNCT__ "MatNorm_SeqAIJ" 1375 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1376 { 1377 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1378 PetscScalar *v = a->a; 1379 PetscReal sum = 0.0; 1380 PetscErrorCode ierr; 1381 PetscInt i,j; 1382 1383 PetscFunctionBegin; 1384 if (type == NORM_FROBENIUS) { 1385 for (i=0; i<a->nz; i++) { 1386 #if defined(PETSC_USE_COMPLEX) 1387 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1388 #else 1389 sum += (*v)*(*v); v++; 1390 #endif 1391 } 1392 *nrm = sqrt(sum); 1393 } else if (type == NORM_1) { 1394 PetscReal *tmp; 1395 PetscInt *jj = a->j; 1396 ierr = PetscMalloc((A->cmap.n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1397 ierr = PetscMemzero(tmp,A->cmap.n*sizeof(PetscReal));CHKERRQ(ierr); 1398 *nrm = 0.0; 1399 for (j=0; j<a->nz; j++) { 1400 tmp[*jj++] += PetscAbsScalar(*v); v++; 1401 } 1402 for (j=0; j<A->cmap.n; j++) { 1403 if (tmp[j] > *nrm) *nrm = tmp[j]; 1404 } 1405 ierr = PetscFree(tmp);CHKERRQ(ierr); 1406 } else if (type == NORM_INFINITY) { 1407 *nrm = 0.0; 1408 for (j=0; j<A->rmap.n; j++) { 1409 v = a->a + a->i[j]; 1410 sum = 0.0; 1411 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1412 sum += PetscAbsScalar(*v); v++; 1413 } 1414 if (sum > *nrm) *nrm = sum; 1415 } 1416 } else { 1417 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1418 } 1419 PetscFunctionReturn(0); 1420 } 1421 1422 #undef __FUNCT__ 1423 #define __FUNCT__ "MatTranspose_SeqAIJ" 1424 PetscErrorCode MatTranspose_SeqAIJ(Mat A,Mat *B) 1425 { 1426 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1427 Mat C; 1428 PetscErrorCode ierr; 1429 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap.n,len,*col; 1430 PetscScalar *array = a->a; 1431 1432 PetscFunctionBegin; 1433 if (!B && m != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1434 ierr = PetscMalloc((1+A->cmap.n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1435 ierr = PetscMemzero(col,(1+A->cmap.n)*sizeof(PetscInt));CHKERRQ(ierr); 1436 1437 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1438 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1439 ierr = MatSetSizes(C,A->cmap.n,m,A->cmap.n,m);CHKERRQ(ierr); 1440 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1441 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1442 ierr = PetscFree(col);CHKERRQ(ierr); 1443 for (i=0; i<m; i++) { 1444 len = ai[i+1]-ai[i]; 1445 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1446 array += len; 1447 aj += len; 1448 } 1449 1450 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1451 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1452 1453 if (B) { 1454 *B = C; 1455 } else { 1456 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1457 } 1458 PetscFunctionReturn(0); 1459 } 1460 1461 EXTERN_C_BEGIN 1462 #undef __FUNCT__ 1463 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1464 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1465 { 1466 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1467 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1468 PetscErrorCode ierr; 1469 PetscInt ma,na,mb,nb, i; 1470 1471 PetscFunctionBegin; 1472 bij = (Mat_SeqAIJ *) B->data; 1473 1474 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1475 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1476 if (ma!=nb || na!=mb){ 1477 *f = PETSC_FALSE; 1478 PetscFunctionReturn(0); 1479 } 1480 aii = aij->i; bii = bij->i; 1481 adx = aij->j; bdx = bij->j; 1482 va = aij->a; vb = bij->a; 1483 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1484 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1485 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1486 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1487 1488 *f = PETSC_TRUE; 1489 for (i=0; i<ma; i++) { 1490 while (aptr[i]<aii[i+1]) { 1491 PetscInt idc,idr; 1492 PetscScalar vc,vr; 1493 /* column/row index/value */ 1494 idc = adx[aptr[i]]; 1495 idr = bdx[bptr[idc]]; 1496 vc = va[aptr[i]]; 1497 vr = vb[bptr[idc]]; 1498 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1499 *f = PETSC_FALSE; 1500 goto done; 1501 } else { 1502 aptr[i]++; 1503 if (B || i!=idc) bptr[idc]++; 1504 } 1505 } 1506 } 1507 done: 1508 ierr = PetscFree(aptr);CHKERRQ(ierr); 1509 if (B) { 1510 ierr = PetscFree(bptr);CHKERRQ(ierr); 1511 } 1512 PetscFunctionReturn(0); 1513 } 1514 EXTERN_C_END 1515 1516 EXTERN_C_BEGIN 1517 #undef __FUNCT__ 1518 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1519 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1520 { 1521 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1522 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; PetscScalar *va,*vb; 1523 PetscErrorCode ierr; 1524 PetscInt ma,na,mb,nb, i; 1525 1526 PetscFunctionBegin; 1527 bij = (Mat_SeqAIJ *) B->data; 1528 1529 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1530 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1531 if (ma!=nb || na!=mb){ 1532 *f = PETSC_FALSE; 1533 PetscFunctionReturn(0); 1534 } 1535 aii = aij->i; bii = bij->i; 1536 adx = aij->j; bdx = bij->j; 1537 va = aij->a; vb = bij->a; 1538 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1539 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1540 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1541 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1542 1543 *f = PETSC_TRUE; 1544 for (i=0; i<ma; i++) { 1545 while (aptr[i]<aii[i+1]) { 1546 PetscInt idc,idr; 1547 PetscScalar vc,vr; 1548 /* column/row index/value */ 1549 idc = adx[aptr[i]]; 1550 idr = bdx[bptr[idc]]; 1551 vc = va[aptr[i]]; 1552 vr = vb[bptr[idc]]; 1553 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1554 *f = PETSC_FALSE; 1555 goto done; 1556 } else { 1557 aptr[i]++; 1558 if (B || i!=idc) bptr[idc]++; 1559 } 1560 } 1561 } 1562 done: 1563 ierr = PetscFree(aptr);CHKERRQ(ierr); 1564 if (B) { 1565 ierr = PetscFree(bptr);CHKERRQ(ierr); 1566 } 1567 PetscFunctionReturn(0); 1568 } 1569 EXTERN_C_END 1570 1571 #undef __FUNCT__ 1572 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1573 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1574 { 1575 PetscErrorCode ierr; 1576 PetscFunctionBegin; 1577 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1578 PetscFunctionReturn(0); 1579 } 1580 1581 #undef __FUNCT__ 1582 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1583 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1584 { 1585 PetscErrorCode ierr; 1586 PetscFunctionBegin; 1587 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1588 PetscFunctionReturn(0); 1589 } 1590 1591 #undef __FUNCT__ 1592 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1593 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1594 { 1595 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1596 PetscScalar *l,*r,x,*v; 1597 PetscErrorCode ierr; 1598 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,M,nz = a->nz,*jj; 1599 1600 PetscFunctionBegin; 1601 if (ll) { 1602 /* The local size is used so that VecMPI can be passed to this routine 1603 by MatDiagonalScale_MPIAIJ */ 1604 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1605 if (m != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1606 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1607 v = a->a; 1608 for (i=0; i<m; i++) { 1609 x = l[i]; 1610 M = a->i[i+1] - a->i[i]; 1611 for (j=0; j<M; j++) { (*v++) *= x;} 1612 } 1613 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1614 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1615 } 1616 if (rr) { 1617 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1618 if (n != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1619 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1620 v = a->a; jj = a->j; 1621 for (i=0; i<nz; i++) { 1622 (*v++) *= r[*jj++]; 1623 } 1624 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1625 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1626 } 1627 PetscFunctionReturn(0); 1628 } 1629 1630 #undef __FUNCT__ 1631 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1632 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1633 { 1634 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1635 PetscErrorCode ierr; 1636 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap.n,*lens; 1637 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1638 PetscInt *irow,*icol,nrows,ncols; 1639 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1640 PetscScalar *a_new,*mat_a; 1641 Mat C; 1642 PetscTruth stride; 1643 1644 PetscFunctionBegin; 1645 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1646 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1647 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1648 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1649 1650 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1651 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1652 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1653 1654 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1655 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1656 if (stride && step == 1) { 1657 /* special case of contiguous rows */ 1658 ierr = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1659 starts = lens + nrows; 1660 /* loop over new rows determining lens and starting points */ 1661 for (i=0; i<nrows; i++) { 1662 kstart = ai[irow[i]]; 1663 kend = kstart + ailen[irow[i]]; 1664 for (k=kstart; k<kend; k++) { 1665 if (aj[k] >= first) { 1666 starts[i] = k; 1667 break; 1668 } 1669 } 1670 sum = 0; 1671 while (k < kend) { 1672 if (aj[k++] >= first+ncols) break; 1673 sum++; 1674 } 1675 lens[i] = sum; 1676 } 1677 /* create submatrix */ 1678 if (scall == MAT_REUSE_MATRIX) { 1679 PetscInt n_cols,n_rows; 1680 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1681 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1682 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1683 C = *B; 1684 } else { 1685 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1686 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1687 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1688 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1689 } 1690 c = (Mat_SeqAIJ*)C->data; 1691 1692 /* loop over rows inserting into submatrix */ 1693 a_new = c->a; 1694 j_new = c->j; 1695 i_new = c->i; 1696 1697 for (i=0; i<nrows; i++) { 1698 ii = starts[i]; 1699 lensi = lens[i]; 1700 for (k=0; k<lensi; k++) { 1701 *j_new++ = aj[ii+k] - first; 1702 } 1703 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1704 a_new += lensi; 1705 i_new[i+1] = i_new[i] + lensi; 1706 c->ilen[i] = lensi; 1707 } 1708 ierr = PetscFree(lens);CHKERRQ(ierr); 1709 } else { 1710 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1711 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1712 1713 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1714 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1715 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1716 /* determine lens of each row */ 1717 for (i=0; i<nrows; i++) { 1718 kstart = ai[irow[i]]; 1719 kend = kstart + a->ilen[irow[i]]; 1720 lens[i] = 0; 1721 for (k=kstart; k<kend; k++) { 1722 if (smap[aj[k]]) { 1723 lens[i]++; 1724 } 1725 } 1726 } 1727 /* Create and fill new matrix */ 1728 if (scall == MAT_REUSE_MATRIX) { 1729 PetscTruth equal; 1730 1731 c = (Mat_SeqAIJ *)((*B)->data); 1732 if ((*B)->rmap.n != nrows || (*B)->cmap.n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1733 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap.n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1734 if (!equal) { 1735 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1736 } 1737 ierr = PetscMemzero(c->ilen,(*B)->rmap.n*sizeof(PetscInt));CHKERRQ(ierr); 1738 C = *B; 1739 } else { 1740 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 1741 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1742 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 1743 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1744 } 1745 c = (Mat_SeqAIJ *)(C->data); 1746 for (i=0; i<nrows; i++) { 1747 row = irow[i]; 1748 kstart = ai[row]; 1749 kend = kstart + a->ilen[row]; 1750 mat_i = c->i[i]; 1751 mat_j = c->j + mat_i; 1752 mat_a = c->a + mat_i; 1753 mat_ilen = c->ilen + i; 1754 for (k=kstart; k<kend; k++) { 1755 if ((tcol=smap[a->j[k]])) { 1756 *mat_j++ = tcol - 1; 1757 *mat_a++ = a->a[k]; 1758 (*mat_ilen)++; 1759 1760 } 1761 } 1762 } 1763 /* Free work space */ 1764 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1765 ierr = PetscFree(smap);CHKERRQ(ierr); 1766 ierr = PetscFree(lens);CHKERRQ(ierr); 1767 } 1768 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1769 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1770 1771 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1772 *B = C; 1773 PetscFunctionReturn(0); 1774 } 1775 1776 /* 1777 */ 1778 #undef __FUNCT__ 1779 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1780 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatFactorInfo *info) 1781 { 1782 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1783 PetscErrorCode ierr; 1784 Mat outA; 1785 PetscTruth row_identity,col_identity; 1786 1787 PetscFunctionBegin; 1788 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1789 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1790 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1791 1792 outA = inA; 1793 inA->factor = FACTOR_LU; 1794 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1795 if (a->row) { ierr = ISDestroy(a->row); CHKERRQ(ierr);} 1796 a->row = row; 1797 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1798 if (a->col) { ierr = ISDestroy(a->col); CHKERRQ(ierr);} 1799 a->col = col; 1800 1801 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1802 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1803 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1804 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1805 1806 if (!a->solve_work) { /* this matrix may have been factored before */ 1807 ierr = PetscMalloc((inA->rmap.n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1808 ierr = PetscLogObjectMemory(inA, (inA->rmap.n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1809 } 1810 1811 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1812 if (row_identity && col_identity) { 1813 ierr = MatLUFactorNumeric_SeqAIJ(inA,info,&outA);CHKERRQ(ierr); 1814 } else { 1815 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(inA,info,&outA);CHKERRQ(ierr); 1816 } 1817 PetscFunctionReturn(0); 1818 } 1819 1820 #include "petscblaslapack.h" 1821 #undef __FUNCT__ 1822 #define __FUNCT__ "MatScale_SeqAIJ" 1823 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1824 { 1825 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1826 PetscBLASInt bnz = (PetscBLASInt)a->nz,one = 1; 1827 PetscScalar oalpha = alpha; 1828 PetscErrorCode ierr; 1829 1830 1831 PetscFunctionBegin; 1832 BLASscal_(&bnz,&oalpha,a->a,&one); 1833 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1834 PetscFunctionReturn(0); 1835 } 1836 1837 #undef __FUNCT__ 1838 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1839 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1840 { 1841 PetscErrorCode ierr; 1842 PetscInt i; 1843 1844 PetscFunctionBegin; 1845 if (scall == MAT_INITIAL_MATRIX) { 1846 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1847 } 1848 1849 for (i=0; i<n; i++) { 1850 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1851 } 1852 PetscFunctionReturn(0); 1853 } 1854 1855 #undef __FUNCT__ 1856 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1857 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1858 { 1859 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1860 PetscErrorCode ierr; 1861 PetscInt row,i,j,k,l,m,n,*idx,*nidx,isz,val; 1862 PetscInt start,end,*ai,*aj; 1863 PetscBT table; 1864 1865 PetscFunctionBegin; 1866 m = A->rmap.n; 1867 ai = a->i; 1868 aj = a->j; 1869 1870 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1871 1872 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1873 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1874 1875 for (i=0; i<is_max; i++) { 1876 /* Initialize the two local arrays */ 1877 isz = 0; 1878 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1879 1880 /* Extract the indices, assume there can be duplicate entries */ 1881 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1882 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1883 1884 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1885 for (j=0; j<n ; ++j){ 1886 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1887 } 1888 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1889 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1890 1891 k = 0; 1892 for (j=0; j<ov; j++){ /* for each overlap */ 1893 n = isz; 1894 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1895 row = nidx[k]; 1896 start = ai[row]; 1897 end = ai[row+1]; 1898 for (l = start; l<end ; l++){ 1899 val = aj[l] ; 1900 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1901 } 1902 } 1903 } 1904 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1905 } 1906 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1907 ierr = PetscFree(nidx);CHKERRQ(ierr); 1908 PetscFunctionReturn(0); 1909 } 1910 1911 /* -------------------------------------------------------------- */ 1912 #undef __FUNCT__ 1913 #define __FUNCT__ "MatPermute_SeqAIJ" 1914 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1915 { 1916 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1917 PetscErrorCode ierr; 1918 PetscInt i,nz,m = A->rmap.n,n = A->cmap.n,*col; 1919 PetscInt *row,*cnew,j,*lens; 1920 IS icolp,irowp; 1921 PetscInt *cwork; 1922 PetscScalar *vwork; 1923 1924 PetscFunctionBegin; 1925 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1926 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1927 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1928 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1929 1930 /* determine lengths of permuted rows */ 1931 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1932 for (i=0; i<m; i++) { 1933 lens[row[i]] = a->i[i+1] - a->i[i]; 1934 } 1935 ierr = MatCreate(A->comm,B);CHKERRQ(ierr); 1936 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 1937 ierr = MatSetType(*B,A->type_name);CHKERRQ(ierr); 1938 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 1939 ierr = PetscFree(lens);CHKERRQ(ierr); 1940 1941 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 1942 for (i=0; i<m; i++) { 1943 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1944 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1945 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1946 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1947 } 1948 ierr = PetscFree(cnew);CHKERRQ(ierr); 1949 (*B)->assembled = PETSC_FALSE; 1950 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1951 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1952 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1953 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1954 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1955 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1956 PetscFunctionReturn(0); 1957 } 1958 1959 #undef __FUNCT__ 1960 #define __FUNCT__ "MatCopy_SeqAIJ" 1961 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1962 { 1963 PetscErrorCode ierr; 1964 1965 PetscFunctionBegin; 1966 /* If the two matrices have the same copy implementation, use fast copy. */ 1967 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1968 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1969 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1970 1971 if (a->i[A->rmap.n] != b->i[B->rmap.n]) { 1972 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1973 } 1974 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap.n])*sizeof(PetscScalar));CHKERRQ(ierr); 1975 } else { 1976 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1977 } 1978 PetscFunctionReturn(0); 1979 } 1980 1981 #undef __FUNCT__ 1982 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1983 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 1984 { 1985 PetscErrorCode ierr; 1986 1987 PetscFunctionBegin; 1988 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1989 PetscFunctionReturn(0); 1990 } 1991 1992 #undef __FUNCT__ 1993 #define __FUNCT__ "MatGetArray_SeqAIJ" 1994 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 1995 { 1996 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1997 PetscFunctionBegin; 1998 *array = a->a; 1999 PetscFunctionReturn(0); 2000 } 2001 2002 #undef __FUNCT__ 2003 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2004 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2005 { 2006 PetscFunctionBegin; 2007 PetscFunctionReturn(0); 2008 } 2009 2010 #undef __FUNCT__ 2011 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2012 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2013 { 2014 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2015 PetscErrorCode ierr; 2016 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2017 PetscScalar dx,*y,*xx,*w3_array; 2018 PetscScalar *vscale_array; 2019 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2020 Vec w1,w2,w3; 2021 void *fctx = coloring->fctx; 2022 PetscTruth flg; 2023 2024 PetscFunctionBegin; 2025 if (!coloring->w1) { 2026 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2027 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2028 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2029 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2030 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2031 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2032 } 2033 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2034 2035 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2036 ierr = PetscOptionsHasName(coloring->prefix,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 2037 if (flg) { 2038 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2039 } else { 2040 PetscTruth assembled; 2041 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2042 if (assembled) { 2043 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2044 } 2045 } 2046 2047 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2048 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2049 2050 /* 2051 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2052 coloring->F for the coarser grids from the finest 2053 */ 2054 if (coloring->F) { 2055 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2056 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2057 if (m1 != m2) { 2058 coloring->F = 0; 2059 } 2060 } 2061 2062 if (coloring->F) { 2063 w1 = coloring->F; 2064 coloring->F = 0; 2065 } else { 2066 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2067 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2068 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2069 } 2070 2071 /* 2072 Compute all the scale factors and share with other processors 2073 */ 2074 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2075 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2076 for (k=0; k<coloring->ncolors; k++) { 2077 /* 2078 Loop over each column associated with color adding the 2079 perturbation to the vector w3. 2080 */ 2081 for (l=0; l<coloring->ncolumns[k]; l++) { 2082 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2083 dx = xx[col]; 2084 if (dx == 0.0) dx = 1.0; 2085 #if !defined(PETSC_USE_COMPLEX) 2086 if (dx < umin && dx >= 0.0) dx = umin; 2087 else if (dx < 0.0 && dx > -umin) dx = -umin; 2088 #else 2089 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2090 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2091 #endif 2092 dx *= epsilon; 2093 vscale_array[col] = 1.0/dx; 2094 } 2095 } 2096 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2097 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2098 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2099 2100 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2101 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2102 2103 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2104 else vscaleforrow = coloring->columnsforrow; 2105 2106 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2107 /* 2108 Loop over each color 2109 */ 2110 for (k=0; k<coloring->ncolors; k++) { 2111 coloring->currentcolor = k; 2112 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2113 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2114 /* 2115 Loop over each column associated with color adding the 2116 perturbation to the vector w3. 2117 */ 2118 for (l=0; l<coloring->ncolumns[k]; l++) { 2119 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2120 dx = xx[col]; 2121 if (dx == 0.0) dx = 1.0; 2122 #if !defined(PETSC_USE_COMPLEX) 2123 if (dx < umin && dx >= 0.0) dx = umin; 2124 else if (dx < 0.0 && dx > -umin) dx = -umin; 2125 #else 2126 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2127 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2128 #endif 2129 dx *= epsilon; 2130 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2131 w3_array[col] += dx; 2132 } 2133 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2134 2135 /* 2136 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2137 */ 2138 2139 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2140 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2141 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2142 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2143 2144 /* 2145 Loop over rows of vector, putting results into Jacobian matrix 2146 */ 2147 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2148 for (l=0; l<coloring->nrows[k]; l++) { 2149 row = coloring->rows[k][l]; 2150 col = coloring->columnsforrow[k][l]; 2151 y[row] *= vscale_array[vscaleforrow[k][l]]; 2152 srow = row + start; 2153 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2154 } 2155 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2156 } 2157 coloring->currentcolor = k; 2158 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2159 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2160 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2161 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2162 PetscFunctionReturn(0); 2163 } 2164 2165 #include "petscblaslapack.h" 2166 #undef __FUNCT__ 2167 #define __FUNCT__ "MatAXPY_SeqAIJ" 2168 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2169 { 2170 PetscErrorCode ierr; 2171 PetscInt i; 2172 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2173 PetscBLASInt one=1,bnz = (PetscBLASInt)x->nz; 2174 2175 PetscFunctionBegin; 2176 if (str == SAME_NONZERO_PATTERN) { 2177 PetscScalar alpha = a; 2178 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2179 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2180 if (y->xtoy && y->XtoY != X) { 2181 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2182 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2183 } 2184 if (!y->xtoy) { /* get xtoy */ 2185 ierr = MatAXPYGetxtoy_Private(X->rmap.n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2186 y->XtoY = X; 2187 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2188 } 2189 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2190 ierr = PetscInfo3(0,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2191 } else { 2192 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2193 } 2194 PetscFunctionReturn(0); 2195 } 2196 2197 #undef __FUNCT__ 2198 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2199 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2200 { 2201 PetscFunctionBegin; 2202 PetscFunctionReturn(0); 2203 } 2204 2205 #undef __FUNCT__ 2206 #define __FUNCT__ "MatConjugate_SeqAIJ" 2207 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2208 { 2209 #if defined(PETSC_USE_COMPLEX) 2210 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2211 PetscInt i,nz; 2212 PetscScalar *a; 2213 2214 PetscFunctionBegin; 2215 nz = aij->nz; 2216 a = aij->a; 2217 for (i=0; i<nz; i++) { 2218 a[i] = PetscConj(a[i]); 2219 } 2220 #else 2221 PetscFunctionBegin; 2222 #endif 2223 PetscFunctionReturn(0); 2224 } 2225 2226 #undef __FUNCT__ 2227 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2228 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2229 { 2230 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2231 PetscErrorCode ierr; 2232 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2233 PetscReal atmp; 2234 PetscScalar *x; 2235 MatScalar *aa; 2236 2237 PetscFunctionBegin; 2238 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2239 aa = a->a; 2240 ai = a->i; 2241 aj = a->j; 2242 2243 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2244 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2245 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2246 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2247 for (i=0; i<m; i++) { 2248 ncols = ai[1] - ai[0]; ai++; 2249 x[i] = 0.0; if (idx) idx[i] = 0; 2250 for (j=0; j<ncols; j++){ 2251 atmp = PetscAbsScalar(*aa); 2252 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2253 aa++; aj++; 2254 } 2255 } 2256 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2257 PetscFunctionReturn(0); 2258 } 2259 2260 #undef __FUNCT__ 2261 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2262 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2263 { 2264 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2265 PetscErrorCode ierr; 2266 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2267 PetscScalar *x; 2268 MatScalar *aa; 2269 2270 PetscFunctionBegin; 2271 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2272 aa = a->a; 2273 ai = a->i; 2274 aj = a->j; 2275 2276 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2277 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2278 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2279 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2280 for (i=0; i<m; i++) { 2281 ncols = ai[1] - ai[0]; ai++; 2282 if (ncols == A->cmap.n) { /* row is dense */ 2283 x[i] = *aa; if (idx) idx[i] = 0; 2284 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2285 x[i] = 0.0; 2286 if (idx) { 2287 idx[i] = 0; /* in case ncols is zero */ 2288 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2289 if (aj[j] > j) { 2290 idx[i] = j; 2291 break; 2292 } 2293 } 2294 } 2295 } 2296 for (j=0; j<ncols; j++){ 2297 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2298 aa++; aj++; 2299 } 2300 } 2301 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2302 PetscFunctionReturn(0); 2303 } 2304 2305 #undef __FUNCT__ 2306 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2307 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2308 { 2309 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2310 PetscErrorCode ierr; 2311 PetscInt i,j,m = A->rmap.n,*ai,*aj,ncols,n; 2312 PetscScalar *x; 2313 MatScalar *aa; 2314 2315 PetscFunctionBegin; 2316 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2317 aa = a->a; 2318 ai = a->i; 2319 aj = a->j; 2320 2321 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2322 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2323 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2324 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2325 for (i=0; i<m; i++) { 2326 ncols = ai[1] - ai[0]; ai++; 2327 if (ncols == A->cmap.n) { /* row is dense */ 2328 x[i] = *aa; if (idx) idx[i] = 0; 2329 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2330 x[i] = 0.0; 2331 if (idx) { /* find first implicit 0.0 in the row */ 2332 idx[i] = 0; /* in case ncols is zero */ 2333 for (j=0;j<ncols;j++) { 2334 if (aj[j] > j) { 2335 idx[i] = j; 2336 break; 2337 } 2338 } 2339 } 2340 } 2341 for (j=0; j<ncols; j++){ 2342 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2343 aa++; aj++; 2344 } 2345 } 2346 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2347 PetscFunctionReturn(0); 2348 } 2349 2350 /* -------------------------------------------------------------------*/ 2351 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2352 MatGetRow_SeqAIJ, 2353 MatRestoreRow_SeqAIJ, 2354 MatMult_SeqAIJ, 2355 /* 4*/ MatMultAdd_SeqAIJ, 2356 MatMultTranspose_SeqAIJ, 2357 MatMultTransposeAdd_SeqAIJ, 2358 MatSolve_SeqAIJ, 2359 MatSolveAdd_SeqAIJ, 2360 MatSolveTranspose_SeqAIJ, 2361 /*10*/ MatSolveTransposeAdd_SeqAIJ, 2362 MatLUFactor_SeqAIJ, 2363 0, 2364 MatRelax_SeqAIJ, 2365 MatTranspose_SeqAIJ, 2366 /*15*/ MatGetInfo_SeqAIJ, 2367 MatEqual_SeqAIJ, 2368 MatGetDiagonal_SeqAIJ, 2369 MatDiagonalScale_SeqAIJ, 2370 MatNorm_SeqAIJ, 2371 /*20*/ 0, 2372 MatAssemblyEnd_SeqAIJ, 2373 MatCompress_SeqAIJ, 2374 MatSetOption_SeqAIJ, 2375 MatZeroEntries_SeqAIJ, 2376 /*25*/ MatZeroRows_SeqAIJ, 2377 MatLUFactorSymbolic_SeqAIJ, 2378 MatLUFactorNumeric_SeqAIJ, 2379 MatCholeskyFactorSymbolic_SeqAIJ, 2380 MatCholeskyFactorNumeric_SeqAIJ, 2381 /*30*/ MatSetUpPreallocation_SeqAIJ, 2382 MatILUFactorSymbolic_SeqAIJ, 2383 MatICCFactorSymbolic_SeqAIJ, 2384 MatGetArray_SeqAIJ, 2385 MatRestoreArray_SeqAIJ, 2386 /*35*/ MatDuplicate_SeqAIJ, 2387 0, 2388 0, 2389 MatILUFactor_SeqAIJ, 2390 0, 2391 /*40*/ MatAXPY_SeqAIJ, 2392 MatGetSubMatrices_SeqAIJ, 2393 MatIncreaseOverlap_SeqAIJ, 2394 MatGetValues_SeqAIJ, 2395 MatCopy_SeqAIJ, 2396 /*45*/ MatGetRowMax_SeqAIJ, 2397 MatScale_SeqAIJ, 2398 0, 2399 MatDiagonalSet_SeqAIJ, 2400 MatILUDTFactor_SeqAIJ, 2401 /*50*/ MatSetBlockSize_SeqAIJ, 2402 MatGetRowIJ_SeqAIJ, 2403 MatRestoreRowIJ_SeqAIJ, 2404 MatGetColumnIJ_SeqAIJ, 2405 MatRestoreColumnIJ_SeqAIJ, 2406 /*55*/ MatFDColoringCreate_SeqAIJ, 2407 0, 2408 0, 2409 MatPermute_SeqAIJ, 2410 0, 2411 /*60*/ 0, 2412 MatDestroy_SeqAIJ, 2413 MatView_SeqAIJ, 2414 0, 2415 0, 2416 /*65*/ 0, 2417 0, 2418 0, 2419 0, 2420 0, 2421 /*70*/ MatGetRowMaxAbs_SeqAIJ, 2422 0, 2423 MatSetColoring_SeqAIJ, 2424 #if defined(PETSC_HAVE_ADIC) 2425 MatSetValuesAdic_SeqAIJ, 2426 #else 2427 0, 2428 #endif 2429 MatSetValuesAdifor_SeqAIJ, 2430 /*75*/ 0, 2431 0, 2432 0, 2433 0, 2434 0, 2435 /*80*/ 0, 2436 0, 2437 0, 2438 0, 2439 MatLoad_SeqAIJ, 2440 /*85*/ MatIsSymmetric_SeqAIJ, 2441 MatIsHermitian_SeqAIJ, 2442 0, 2443 0, 2444 0, 2445 /*90*/ MatMatMult_SeqAIJ_SeqAIJ, 2446 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2447 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2448 MatPtAP_Basic, 2449 MatPtAPSymbolic_SeqAIJ, 2450 /*95*/ MatPtAPNumeric_SeqAIJ, 2451 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2452 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2453 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2454 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2455 /*100*/MatPtAPNumeric_SeqAIJ_SeqAIJ, 2456 0, 2457 0, 2458 MatConjugate_SeqAIJ, 2459 0, 2460 /*105*/MatSetValuesRow_SeqAIJ, 2461 MatRealPart_SeqAIJ, 2462 MatImaginaryPart_SeqAIJ, 2463 0, 2464 0, 2465 /*110*/MatMatSolve_SeqAIJ, 2466 0, 2467 MatGetRowMin_SeqAIJ 2468 }; 2469 2470 EXTERN_C_BEGIN 2471 #undef __FUNCT__ 2472 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2473 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2474 { 2475 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2476 PetscInt i,nz,n; 2477 2478 PetscFunctionBegin; 2479 2480 nz = aij->maxnz; 2481 n = mat->rmap.n; 2482 for (i=0; i<nz; i++) { 2483 aij->j[i] = indices[i]; 2484 } 2485 aij->nz = nz; 2486 for (i=0; i<n; i++) { 2487 aij->ilen[i] = aij->imax[i]; 2488 } 2489 2490 PetscFunctionReturn(0); 2491 } 2492 EXTERN_C_END 2493 2494 #undef __FUNCT__ 2495 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2496 /*@ 2497 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2498 in the matrix. 2499 2500 Input Parameters: 2501 + mat - the SeqAIJ matrix 2502 - indices - the column indices 2503 2504 Level: advanced 2505 2506 Notes: 2507 This can be called if you have precomputed the nonzero structure of the 2508 matrix and want to provide it to the matrix object to improve the performance 2509 of the MatSetValues() operation. 2510 2511 You MUST have set the correct numbers of nonzeros per row in the call to 2512 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2513 2514 MUST be called before any calls to MatSetValues(); 2515 2516 The indices should start with zero, not one. 2517 2518 @*/ 2519 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2520 { 2521 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2522 2523 PetscFunctionBegin; 2524 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2525 PetscValidPointer(indices,2); 2526 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2527 if (f) { 2528 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2529 } else { 2530 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2531 } 2532 PetscFunctionReturn(0); 2533 } 2534 2535 /* ----------------------------------------------------------------------------------------*/ 2536 2537 EXTERN_C_BEGIN 2538 #undef __FUNCT__ 2539 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2540 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2541 { 2542 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2543 PetscErrorCode ierr; 2544 size_t nz = aij->i[mat->rmap.n]; 2545 2546 PetscFunctionBegin; 2547 if (aij->nonew != 1) { 2548 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2549 } 2550 2551 /* allocate space for values if not already there */ 2552 if (!aij->saved_values) { 2553 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2554 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2555 } 2556 2557 /* copy values over */ 2558 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2559 PetscFunctionReturn(0); 2560 } 2561 EXTERN_C_END 2562 2563 #undef __FUNCT__ 2564 #define __FUNCT__ "MatStoreValues" 2565 /*@ 2566 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2567 example, reuse of the linear part of a Jacobian, while recomputing the 2568 nonlinear portion. 2569 2570 Collect on Mat 2571 2572 Input Parameters: 2573 . mat - the matrix (currently only AIJ matrices support this option) 2574 2575 Level: advanced 2576 2577 Common Usage, with SNESSolve(): 2578 $ Create Jacobian matrix 2579 $ Set linear terms into matrix 2580 $ Apply boundary conditions to matrix, at this time matrix must have 2581 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2582 $ boundary conditions again will not change the nonzero structure 2583 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2584 $ ierr = MatStoreValues(mat); 2585 $ Call SNESSetJacobian() with matrix 2586 $ In your Jacobian routine 2587 $ ierr = MatRetrieveValues(mat); 2588 $ Set nonlinear terms in matrix 2589 2590 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2591 $ // build linear portion of Jacobian 2592 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2593 $ ierr = MatStoreValues(mat); 2594 $ loop over nonlinear iterations 2595 $ ierr = MatRetrieveValues(mat); 2596 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2597 $ // call MatAssemblyBegin/End() on matrix 2598 $ Solve linear system with Jacobian 2599 $ endloop 2600 2601 Notes: 2602 Matrix must already be assemblied before calling this routine 2603 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2604 calling this routine. 2605 2606 When this is called multiple times it overwrites the previous set of stored values 2607 and does not allocated additional space. 2608 2609 .seealso: MatRetrieveValues() 2610 2611 @*/ 2612 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2613 { 2614 PetscErrorCode ierr,(*f)(Mat); 2615 2616 PetscFunctionBegin; 2617 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2618 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2619 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2620 2621 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2622 if (f) { 2623 ierr = (*f)(mat);CHKERRQ(ierr); 2624 } else { 2625 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2626 } 2627 PetscFunctionReturn(0); 2628 } 2629 2630 EXTERN_C_BEGIN 2631 #undef __FUNCT__ 2632 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2633 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2634 { 2635 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2636 PetscErrorCode ierr; 2637 PetscInt nz = aij->i[mat->rmap.n]; 2638 2639 PetscFunctionBegin; 2640 if (aij->nonew != 1) { 2641 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2642 } 2643 if (!aij->saved_values) { 2644 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2645 } 2646 /* copy values over */ 2647 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2648 PetscFunctionReturn(0); 2649 } 2650 EXTERN_C_END 2651 2652 #undef __FUNCT__ 2653 #define __FUNCT__ "MatRetrieveValues" 2654 /*@ 2655 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2656 example, reuse of the linear part of a Jacobian, while recomputing the 2657 nonlinear portion. 2658 2659 Collect on Mat 2660 2661 Input Parameters: 2662 . mat - the matrix (currently on AIJ matrices support this option) 2663 2664 Level: advanced 2665 2666 .seealso: MatStoreValues() 2667 2668 @*/ 2669 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2670 { 2671 PetscErrorCode ierr,(*f)(Mat); 2672 2673 PetscFunctionBegin; 2674 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2675 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2676 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2677 2678 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2679 if (f) { 2680 ierr = (*f)(mat);CHKERRQ(ierr); 2681 } else { 2682 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2683 } 2684 PetscFunctionReturn(0); 2685 } 2686 2687 2688 /* --------------------------------------------------------------------------------*/ 2689 #undef __FUNCT__ 2690 #define __FUNCT__ "MatCreateSeqAIJ" 2691 /*@C 2692 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2693 (the default parallel PETSc format). For good matrix assembly performance 2694 the user should preallocate the matrix storage by setting the parameter nz 2695 (or the array nnz). By setting these parameters accurately, performance 2696 during matrix assembly can be increased by more than a factor of 50. 2697 2698 Collective on MPI_Comm 2699 2700 Input Parameters: 2701 + comm - MPI communicator, set to PETSC_COMM_SELF 2702 . m - number of rows 2703 . n - number of columns 2704 . nz - number of nonzeros per row (same for all rows) 2705 - nnz - array containing the number of nonzeros in the various rows 2706 (possibly different for each row) or PETSC_NULL 2707 2708 Output Parameter: 2709 . A - the matrix 2710 2711 Notes: 2712 If nnz is given then nz is ignored 2713 2714 The AIJ format (also called the Yale sparse matrix format or 2715 compressed row storage), is fully compatible with standard Fortran 77 2716 storage. That is, the stored row and column indices can begin at 2717 either one (as in Fortran) or zero. See the users' manual for details. 2718 2719 Specify the preallocated storage with either nz or nnz (not both). 2720 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2721 allocation. For large problems you MUST preallocate memory or you 2722 will get TERRIBLE performance, see the users' manual chapter on matrices. 2723 2724 By default, this format uses inodes (identical nodes) when possible, to 2725 improve numerical efficiency of matrix-vector products and solves. We 2726 search for consecutive rows with the same nonzero structure, thereby 2727 reusing matrix information to achieve increased efficiency. 2728 2729 Options Database Keys: 2730 + -mat_no_inode - Do not use inodes 2731 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2732 - -mat_aij_oneindex - Internally use indexing starting at 1 2733 rather than 0. Note that when calling MatSetValues(), 2734 the user still MUST index entries starting at 0! 2735 2736 Level: intermediate 2737 2738 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2739 2740 @*/ 2741 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2742 { 2743 PetscErrorCode ierr; 2744 2745 PetscFunctionBegin; 2746 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2747 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2748 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2749 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2750 PetscFunctionReturn(0); 2751 } 2752 2753 #undef __FUNCT__ 2754 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2755 /*@C 2756 MatSeqAIJSetPreallocation - For good matrix assembly performance 2757 the user should preallocate the matrix storage by setting the parameter nz 2758 (or the array nnz). By setting these parameters accurately, performance 2759 during matrix assembly can be increased by more than a factor of 50. 2760 2761 Collective on MPI_Comm 2762 2763 Input Parameters: 2764 + B - The matrix 2765 . nz - number of nonzeros per row (same for all rows) 2766 - nnz - array containing the number of nonzeros in the various rows 2767 (possibly different for each row) or PETSC_NULL 2768 2769 Notes: 2770 If nnz is given then nz is ignored 2771 2772 The AIJ format (also called the Yale sparse matrix format or 2773 compressed row storage), is fully compatible with standard Fortran 77 2774 storage. That is, the stored row and column indices can begin at 2775 either one (as in Fortran) or zero. See the users' manual for details. 2776 2777 Specify the preallocated storage with either nz or nnz (not both). 2778 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2779 allocation. For large problems you MUST preallocate memory or you 2780 will get TERRIBLE performance, see the users' manual chapter on matrices. 2781 2782 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2783 entries or columns indices 2784 2785 By default, this format uses inodes (identical nodes) when possible, to 2786 improve numerical efficiency of matrix-vector products and solves. We 2787 search for consecutive rows with the same nonzero structure, thereby 2788 reusing matrix information to achieve increased efficiency. 2789 2790 Options Database Keys: 2791 + -mat_no_inode - Do not use inodes 2792 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2793 - -mat_aij_oneindex - Internally use indexing starting at 1 2794 rather than 0. Note that when calling MatSetValues(), 2795 the user still MUST index entries starting at 0! 2796 2797 Level: intermediate 2798 2799 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2800 2801 @*/ 2802 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2803 { 2804 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2805 2806 PetscFunctionBegin; 2807 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2808 if (f) { 2809 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2810 } 2811 PetscFunctionReturn(0); 2812 } 2813 2814 EXTERN_C_BEGIN 2815 #undef __FUNCT__ 2816 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2817 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2818 { 2819 Mat_SeqAIJ *b; 2820 PetscTruth skipallocation = PETSC_FALSE; 2821 PetscErrorCode ierr; 2822 PetscInt i; 2823 2824 PetscFunctionBegin; 2825 2826 if (nz == MAT_SKIP_ALLOCATION) { 2827 skipallocation = PETSC_TRUE; 2828 nz = 0; 2829 } 2830 2831 B->rmap.bs = B->cmap.bs = 1; 2832 ierr = PetscMapSetUp(&B->rmap);CHKERRQ(ierr); 2833 ierr = PetscMapSetUp(&B->cmap);CHKERRQ(ierr); 2834 2835 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2836 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2837 if (nnz) { 2838 for (i=0; i<B->rmap.n; i++) { 2839 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2840 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); 2841 } 2842 } 2843 2844 B->preallocated = PETSC_TRUE; 2845 b = (Mat_SeqAIJ*)B->data; 2846 2847 if (!skipallocation) { 2848 ierr = PetscMalloc2(B->rmap.n,PetscInt,&b->imax,B->rmap.n,PetscInt,&b->ilen);CHKERRQ(ierr); 2849 ierr = PetscLogObjectMemory(B,2*B->rmap.n*sizeof(PetscInt));CHKERRQ(ierr); 2850 if (!nnz) { 2851 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2852 else if (nz <= 0) nz = 1; 2853 for (i=0; i<B->rmap.n; i++) b->imax[i] = nz; 2854 nz = nz*B->rmap.n; 2855 } else { 2856 nz = 0; 2857 for (i=0; i<B->rmap.n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2858 } 2859 2860 /* b->ilen will count nonzeros in each row so far. */ 2861 for (i=0; i<B->rmap.n; i++) { b->ilen[i] = 0;} 2862 2863 /* allocate the matrix space */ 2864 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap.n+1,PetscInt,&b->i);CHKERRQ(ierr); 2865 ierr = PetscLogObjectMemory(B,(B->rmap.n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 2866 b->i[0] = 0; 2867 for (i=1; i<B->rmap.n+1; i++) { 2868 b->i[i] = b->i[i-1] + b->imax[i-1]; 2869 } 2870 b->singlemalloc = PETSC_TRUE; 2871 b->free_a = PETSC_TRUE; 2872 b->free_ij = PETSC_TRUE; 2873 } else { 2874 b->free_a = PETSC_FALSE; 2875 b->free_ij = PETSC_FALSE; 2876 } 2877 2878 b->nz = 0; 2879 b->maxnz = nz; 2880 B->info.nz_unneeded = (double)b->maxnz; 2881 PetscFunctionReturn(0); 2882 } 2883 EXTERN_C_END 2884 2885 #undef __FUNCT__ 2886 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 2887 /*@C 2888 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 2889 2890 Input Parameters: 2891 + B - the matrix 2892 . i - the indices into j for the start of each row (starts with zero) 2893 . j - the column indices for each row (starts with zero) these must be sorted for each row 2894 - v - optional values in the matrix 2895 2896 Contributed by: Lisandro Dalchin 2897 2898 Level: developer 2899 2900 .keywords: matrix, aij, compressed row, sparse, sequential 2901 2902 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 2903 @*/ 2904 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 2905 { 2906 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 2907 PetscErrorCode ierr; 2908 2909 PetscFunctionBegin; 2910 PetscValidHeaderSpecific(B,MAT_COOKIE,1); 2911 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 2912 if (f) { 2913 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 2914 } 2915 PetscFunctionReturn(0); 2916 } 2917 2918 EXTERN_C_BEGIN 2919 #undef __FUNCT__ 2920 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 2921 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 2922 { 2923 PetscInt i; 2924 PetscInt m,n; 2925 PetscInt nz; 2926 PetscInt *nnz, nz_max = 0; 2927 PetscScalar *values; 2928 PetscErrorCode ierr; 2929 2930 PetscFunctionBegin; 2931 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 2932 2933 if (Ii[0]) { 2934 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 2935 } 2936 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 2937 for(i = 0; i < m; i++) { 2938 nz = Ii[i+1]- Ii[i]; 2939 nz_max = PetscMax(nz_max, nz); 2940 if (nz < 0) { 2941 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 2942 } 2943 nnz[i] = nz; 2944 } 2945 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 2946 ierr = PetscFree(nnz);CHKERRQ(ierr); 2947 2948 if (v) { 2949 values = (PetscScalar*) v; 2950 } else { 2951 ierr = PetscMalloc((nz_max+1)*sizeof(PetscScalar), &values);CHKERRQ(ierr); 2952 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 2953 } 2954 2955 for(i = 0; i < m; i++) { 2956 nz = Ii[i+1] - Ii[i]; 2957 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 2958 } 2959 2960 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2961 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2962 2963 if (!v) { 2964 ierr = PetscFree(values);CHKERRQ(ierr); 2965 } 2966 PetscFunctionReturn(0); 2967 } 2968 EXTERN_C_END 2969 2970 /*MC 2971 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2972 based on compressed sparse row format. 2973 2974 Options Database Keys: 2975 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2976 2977 Level: beginner 2978 2979 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2980 M*/ 2981 2982 EXTERN_C_BEGIN 2983 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 2984 EXTERN_C_END 2985 2986 EXTERN_C_BEGIN 2987 #undef __FUNCT__ 2988 #define __FUNCT__ "MatCreate_SeqAIJ" 2989 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 2990 { 2991 Mat_SeqAIJ *b; 2992 PetscErrorCode ierr; 2993 PetscMPIInt size; 2994 2995 PetscFunctionBegin; 2996 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2997 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2998 2999 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3000 B->data = (void*)b; 3001 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3002 B->factor = 0; 3003 B->mapping = 0; 3004 b->row = 0; 3005 b->col = 0; 3006 b->icol = 0; 3007 b->reallocs = 0; 3008 b->ignorezeroentries = PETSC_FALSE; 3009 b->roworiented = PETSC_TRUE; 3010 b->nonew = 0; 3011 b->diag = 0; 3012 b->solve_work = 0; 3013 B->spptr = 0; 3014 b->saved_values = 0; 3015 b->idiag = 0; 3016 b->ssor = 0; 3017 b->keepzeroedrows = PETSC_FALSE; 3018 b->xtoy = 0; 3019 b->XtoY = 0; 3020 b->compressedrow.use = PETSC_FALSE; 3021 b->compressedrow.nrows = B->rmap.n; 3022 b->compressedrow.i = PETSC_NULL; 3023 b->compressedrow.rindex = PETSC_NULL; 3024 b->compressedrow.checked = PETSC_FALSE; 3025 B->same_nonzero = PETSC_FALSE; 3026 3027 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3028 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3029 "MatSeqAIJSetColumnIndices_SeqAIJ", 3030 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3031 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3032 "MatStoreValues_SeqAIJ", 3033 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3034 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3035 "MatRetrieveValues_SeqAIJ", 3036 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3037 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3038 "MatConvert_SeqAIJ_SeqSBAIJ", 3039 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3040 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3041 "MatConvert_SeqAIJ_SeqBAIJ", 3042 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3043 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3044 "MatConvert_SeqAIJ_SeqCSRPERM", 3045 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3046 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3047 "MatConvert_SeqAIJ_SeqCRL", 3048 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3049 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3050 "MatIsTranspose_SeqAIJ", 3051 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3052 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3053 "MatIsHermitianTranspose_SeqAIJ", 3054 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3055 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3056 "MatSeqAIJSetPreallocation_SeqAIJ", 3057 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3058 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3059 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3060 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3061 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3062 "MatReorderForNonzeroDiagonal_SeqAIJ", 3063 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3064 ierr = MatCreate_Inode(B);CHKERRQ(ierr); 3065 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3066 PetscFunctionReturn(0); 3067 } 3068 EXTERN_C_END 3069 3070 #undef __FUNCT__ 3071 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3072 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3073 { 3074 Mat C; 3075 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3076 PetscErrorCode ierr; 3077 PetscInt i,m = A->rmap.n; 3078 3079 PetscFunctionBegin; 3080 *B = 0; 3081 ierr = MatCreate(A->comm,&C);CHKERRQ(ierr); 3082 ierr = MatSetSizes(C,A->rmap.n,A->cmap.n,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 3083 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 3084 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3085 3086 ierr = PetscMapCopy(A->comm,&A->rmap,&C->rmap);CHKERRQ(ierr); 3087 ierr = PetscMapCopy(A->comm,&A->cmap,&C->cmap);CHKERRQ(ierr); 3088 3089 c = (Mat_SeqAIJ*)C->data; 3090 3091 C->factor = A->factor; 3092 3093 c->row = 0; 3094 c->col = 0; 3095 c->icol = 0; 3096 c->reallocs = 0; 3097 3098 C->assembled = PETSC_TRUE; 3099 3100 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3101 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3102 for (i=0; i<m; i++) { 3103 c->imax[i] = a->imax[i]; 3104 c->ilen[i] = a->ilen[i]; 3105 } 3106 3107 /* allocate the matrix space */ 3108 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3109 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3110 c->singlemalloc = PETSC_TRUE; 3111 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3112 if (m > 0) { 3113 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3114 if (cpvalues == MAT_COPY_VALUES) { 3115 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3116 } else { 3117 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3118 } 3119 } 3120 3121 c->ignorezeroentries = a->ignorezeroentries; 3122 c->roworiented = a->roworiented; 3123 c->nonew = a->nonew; 3124 if (a->diag) { 3125 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3126 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3127 for (i=0; i<m; i++) { 3128 c->diag[i] = a->diag[i]; 3129 } 3130 } else c->diag = 0; 3131 c->solve_work = 0; 3132 c->saved_values = 0; 3133 c->idiag = 0; 3134 c->ssor = 0; 3135 c->keepzeroedrows = a->keepzeroedrows; 3136 c->free_a = PETSC_TRUE; 3137 c->free_ij = PETSC_TRUE; 3138 c->xtoy = 0; 3139 c->XtoY = 0; 3140 3141 c->nz = a->nz; 3142 c->maxnz = a->maxnz; 3143 C->preallocated = PETSC_TRUE; 3144 3145 c->compressedrow.use = a->compressedrow.use; 3146 c->compressedrow.nrows = a->compressedrow.nrows; 3147 c->compressedrow.checked = a->compressedrow.checked; 3148 if ( a->compressedrow.checked && a->compressedrow.use){ 3149 i = a->compressedrow.nrows; 3150 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 3151 c->compressedrow.rindex = c->compressedrow.i + i + 1; 3152 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3153 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3154 } else { 3155 c->compressedrow.use = PETSC_FALSE; 3156 c->compressedrow.i = PETSC_NULL; 3157 c->compressedrow.rindex = PETSC_NULL; 3158 } 3159 C->same_nonzero = A->same_nonzero; 3160 ierr = MatDuplicate_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3161 3162 *B = C; 3163 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 3164 PetscFunctionReturn(0); 3165 } 3166 3167 #undef __FUNCT__ 3168 #define __FUNCT__ "MatLoad_SeqAIJ" 3169 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, MatType type,Mat *A) 3170 { 3171 Mat_SeqAIJ *a; 3172 Mat B; 3173 PetscErrorCode ierr; 3174 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3175 int fd; 3176 PetscMPIInt size; 3177 MPI_Comm comm; 3178 3179 PetscFunctionBegin; 3180 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3181 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3182 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 3183 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3184 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3185 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3186 M = header[1]; N = header[2]; nz = header[3]; 3187 3188 if (nz < 0) { 3189 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3190 } 3191 3192 /* read in row lengths */ 3193 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3194 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3195 3196 /* check if sum of rowlengths is same as nz */ 3197 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3198 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3199 3200 /* create our matrix */ 3201 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3202 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3203 ierr = MatSetType(B,type);CHKERRQ(ierr); 3204 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3205 a = (Mat_SeqAIJ*)B->data; 3206 3207 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3208 3209 /* read in nonzero values */ 3210 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3211 3212 /* set matrix "i" values */ 3213 a->i[0] = 0; 3214 for (i=1; i<= M; i++) { 3215 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3216 a->ilen[i-1] = rowlengths[i-1]; 3217 } 3218 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3219 3220 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3221 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3222 *A = B; 3223 PetscFunctionReturn(0); 3224 } 3225 3226 #undef __FUNCT__ 3227 #define __FUNCT__ "MatEqual_SeqAIJ" 3228 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3229 { 3230 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3231 PetscErrorCode ierr; 3232 3233 PetscFunctionBegin; 3234 /* If the matrix dimensions are not equal,or no of nonzeros */ 3235 if ((A->rmap.n != B->rmap.n) || (A->cmap.n != B->cmap.n) ||(a->nz != b->nz)) { 3236 *flg = PETSC_FALSE; 3237 PetscFunctionReturn(0); 3238 } 3239 3240 /* if the a->i are the same */ 3241 ierr = PetscMemcmp(a->i,b->i,(A->rmap.n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3242 if (!*flg) PetscFunctionReturn(0); 3243 3244 /* if a->j are the same */ 3245 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3246 if (!*flg) PetscFunctionReturn(0); 3247 3248 /* if a->a are the same */ 3249 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3250 3251 PetscFunctionReturn(0); 3252 3253 } 3254 3255 #undef __FUNCT__ 3256 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3257 /*@ 3258 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3259 provided by the user. 3260 3261 Collective on MPI_Comm 3262 3263 Input Parameters: 3264 + comm - must be an MPI communicator of size 1 3265 . m - number of rows 3266 . n - number of columns 3267 . i - row indices 3268 . j - column indices 3269 - a - matrix values 3270 3271 Output Parameter: 3272 . mat - the matrix 3273 3274 Level: intermediate 3275 3276 Notes: 3277 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3278 once the matrix is destroyed 3279 3280 You cannot set new nonzero locations into this matrix, that will generate an error. 3281 3282 The i and j indices are 0 based 3283 3284 The format which is used for the sparse matrix input, is equivalent to a 3285 row-major ordering.. i.e for the following matrix, the input data expected is 3286 as shown: 3287 3288 1 0 0 3289 2 0 3 3290 4 5 6 3291 3292 i = {0,1,3,6} [size = nrow+1 = 3+1] 3293 j = {0,0,2,0,1,2} [size = nz = 6] 3294 v = {1,2,3,4,5,6} [size = nz = 6] 3295 3296 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3297 3298 @*/ 3299 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3300 { 3301 PetscErrorCode ierr; 3302 PetscInt ii; 3303 Mat_SeqAIJ *aij; 3304 3305 PetscFunctionBegin; 3306 if (i[0]) { 3307 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3308 } 3309 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3310 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3311 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3312 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3313 aij = (Mat_SeqAIJ*)(*mat)->data; 3314 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3315 3316 aij->i = i; 3317 aij->j = j; 3318 aij->a = a; 3319 aij->singlemalloc = PETSC_FALSE; 3320 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3321 aij->free_a = PETSC_FALSE; 3322 aij->free_ij = PETSC_FALSE; 3323 3324 for (ii=0; ii<m; ii++) { 3325 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3326 #if defined(PETSC_USE_DEBUG) 3327 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]); 3328 #endif 3329 } 3330 #if defined(PETSC_USE_DEBUG) 3331 for (ii=0; ii<aij->i[m]; ii++) { 3332 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3333 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3334 } 3335 #endif 3336 3337 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3338 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3339 PetscFunctionReturn(0); 3340 } 3341 3342 #undef __FUNCT__ 3343 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3344 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3345 { 3346 PetscErrorCode ierr; 3347 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3348 3349 PetscFunctionBegin; 3350 if (coloring->ctype == IS_COLORING_GLOBAL) { 3351 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3352 a->coloring = coloring; 3353 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3354 PetscInt i,*larray; 3355 ISColoring ocoloring; 3356 ISColoringValue *colors; 3357 3358 /* set coloring for diagonal portion */ 3359 ierr = PetscMalloc((A->cmap.n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3360 for (i=0; i<A->cmap.n; i++) { 3361 larray[i] = i; 3362 } 3363 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap.n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3364 ierr = PetscMalloc((A->cmap.n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3365 for (i=0; i<A->cmap.n; i++) { 3366 colors[i] = coloring->colors[larray[i]]; 3367 } 3368 ierr = PetscFree(larray);CHKERRQ(ierr); 3369 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap.n,colors,&ocoloring);CHKERRQ(ierr); 3370 a->coloring = ocoloring; 3371 } 3372 PetscFunctionReturn(0); 3373 } 3374 3375 #if defined(PETSC_HAVE_ADIC) 3376 EXTERN_C_BEGIN 3377 #include "adic/ad_utils.h" 3378 EXTERN_C_END 3379 3380 #undef __FUNCT__ 3381 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3382 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3383 { 3384 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3385 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3386 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3387 ISColoringValue *color; 3388 3389 PetscFunctionBegin; 3390 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3391 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3392 color = a->coloring->colors; 3393 /* loop over rows */ 3394 for (i=0; i<m; i++) { 3395 nz = ii[i+1] - ii[i]; 3396 /* loop over columns putting computed value into matrix */ 3397 for (j=0; j<nz; j++) { 3398 *v++ = values[color[*jj++]]; 3399 } 3400 values += nlen; /* jump to next row of derivatives */ 3401 } 3402 PetscFunctionReturn(0); 3403 } 3404 #endif 3405 3406 #undef __FUNCT__ 3407 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3408 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3409 { 3410 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3411 PetscInt m = A->rmap.n,*ii = a->i,*jj = a->j,nz,i,j; 3412 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3413 ISColoringValue *color; 3414 3415 PetscFunctionBegin; 3416 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3417 color = a->coloring->colors; 3418 /* loop over rows */ 3419 for (i=0; i<m; i++) { 3420 nz = ii[i+1] - ii[i]; 3421 /* loop over columns putting computed value into matrix */ 3422 for (j=0; j<nz; j++) { 3423 *v++ = values[color[*jj++]]; 3424 } 3425 values += nl; /* jump to next row of derivatives */ 3426 } 3427 PetscFunctionReturn(0); 3428 } 3429 3430 /* 3431 Special version for direct calls from Fortran 3432 */ 3433 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3434 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3435 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3436 #define matsetvaluesseqaij_ matsetvaluesseqaij 3437 #endif 3438 3439 /* Change these macros so can be used in void function */ 3440 #undef CHKERRQ 3441 #define CHKERRQ(ierr) CHKERRABORT(A->comm,ierr) 3442 #undef SETERRQ2 3443 #define SETERRQ2(ierr,b,c,d) CHKERRABORT(A->comm,ierr) 3444 3445 EXTERN_C_BEGIN 3446 #undef __FUNCT__ 3447 #define __FUNCT__ "matsetvaluesseqaij_" 3448 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3449 { 3450 Mat A = *AA; 3451 PetscInt m = *mm, n = *nn; 3452 InsertMode is = *isis; 3453 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3454 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3455 PetscInt *imax,*ai,*ailen; 3456 PetscErrorCode ierr; 3457 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3458 PetscScalar *ap,value,*aa; 3459 PetscTruth ignorezeroentries = a->ignorezeroentries; 3460 PetscTruth roworiented = a->roworiented; 3461 3462 PetscFunctionBegin; 3463 MatPreallocated(A); 3464 imax = a->imax; 3465 ai = a->i; 3466 ailen = a->ilen; 3467 aj = a->j; 3468 aa = a->a; 3469 3470 for (k=0; k<m; k++) { /* loop over added rows */ 3471 row = im[k]; 3472 if (row < 0) continue; 3473 #if defined(PETSC_USE_DEBUG) 3474 if (row >= A->rmap.n) SETERRABORT(A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3475 #endif 3476 rp = aj + ai[row]; ap = aa + ai[row]; 3477 rmax = imax[row]; nrow = ailen[row]; 3478 low = 0; 3479 high = nrow; 3480 for (l=0; l<n; l++) { /* loop over added columns */ 3481 if (in[l] < 0) continue; 3482 #if defined(PETSC_USE_DEBUG) 3483 if (in[l] >= A->cmap.n) SETERRABORT(A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3484 #endif 3485 col = in[l]; 3486 if (roworiented) { 3487 value = v[l + k*n]; 3488 } else { 3489 value = v[k + l*m]; 3490 } 3491 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3492 3493 if (col <= lastcol) low = 0; else high = nrow; 3494 lastcol = col; 3495 while (high-low > 5) { 3496 t = (low+high)/2; 3497 if (rp[t] > col) high = t; 3498 else low = t; 3499 } 3500 for (i=low; i<high; i++) { 3501 if (rp[i] > col) break; 3502 if (rp[i] == col) { 3503 if (is == ADD_VALUES) ap[i] += value; 3504 else ap[i] = value; 3505 goto noinsert; 3506 } 3507 } 3508 if (value == 0.0 && ignorezeroentries) goto noinsert; 3509 if (nonew == 1) goto noinsert; 3510 if (nonew == -1) SETERRABORT(A->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3511 MatSeqXAIJReallocateAIJ(A,A->rmap.n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3512 N = nrow++ - 1; a->nz++; high++; 3513 /* shift up all the later entries in this row */ 3514 for (ii=N; ii>=i; ii--) { 3515 rp[ii+1] = rp[ii]; 3516 ap[ii+1] = ap[ii]; 3517 } 3518 rp[i] = col; 3519 ap[i] = value; 3520 noinsert:; 3521 low = i + 1; 3522 } 3523 ailen[row] = nrow; 3524 } 3525 A->same_nonzero = PETSC_FALSE; 3526 PetscFunctionReturnVoid(); 3527 } 3528 EXTERN_C_END 3529