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