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