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