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