| #define __FUNCT__ "MatMult_Linear_H1A" |
| PetscErrorCode MatMult_Linear_H1A(Mat A,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_H1B" |
| PetscErrorCode MatMult_Linear_H1B(Mat B,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(B,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_H2A" |
| PetscErrorCode MatMult_Linear_H2A(Mat A,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_H2B" |
| PetscErrorCode MatMult_Linear_H2B(Mat B,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(B,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_N1A" |
| PetscErrorCode MatMult_Linear_N1A(Mat A,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_N1B" |
| PetscErrorCode MatMult_Linear_N1B(Mat B,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(B,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_N2A" |
| PetscErrorCode MatMult_Linear_N2A(Mat A,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_N2B" |
| PetscErrorCode MatMult_Linear_N2B(Mat B,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(B,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_S1A" |
| PetscErrorCode MatMult_Linear_S1A(Mat A,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_S1B" |
| PetscErrorCode MatMult_Linear_S1B(Mat B,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(B,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_S2A" |
| PetscErrorCode MatMult_Linear_S2A(Mat A,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatMult_Linear_S2B" |
| PetscErrorCode MatMult_Linear_S2B(Mat B,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| QEP_LINEAR *ctx; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(B,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->M,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| @*/ |
| PetscErrorCode STComputeExplicitOperator(ST st,Mat *mat) |
| { |
| PetscErrorCode ierr; |
| Vec in,out; |
| PetscInt i,M,m,*rows,start,end; |
| PetscScalar *array,one = 1.0; |
| PetscErrorCode ierr; |
| Vec in,out; |
| PetscInt i,M,m,*rows,start,end; |
| const PetscScalar *array; |
| PetscScalar one = 1.0; |
| PetscFunctionBegin; |
| PetscValidHeaderSpecific(st,ST_CLASSID,1); |
| ierr = STApply(st,in,out);CHKERRQ(ierr); |
| ierr = VecGetArray(out,&array);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(out,&array);CHKERRQ(ierr); |
| ierr = MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); |
| ierr = VecRestoreArray(out,&array);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(out,&array);CHKERRQ(ierr); |
| } |
| ierr = PetscFree(rows);CHKERRQ(ierr); |
| ierr = VecDestroy(&in);CHKERRQ(ierr); |
| */ |
| PetscErrorCode EPSBuildBalance_Krylov(EPS eps) |
| { |
| Vec z,p,r; |
| PetscInt i,j; |
| PetscReal norma; |
| PetscScalar *pz,*pr,*pp,*pD; |
| PetscErrorCode ierr; |
| Vec z,p,r; |
| PetscInt i,j; |
| PetscReal norma; |
| PetscScalar *pz,*pD; |
| const PetscScalar *pr,*pp; |
| PetscErrorCode ierr; |
| PetscFunctionBegin; |
| ierr = VecDuplicate(eps->V[0],&r);CHKERRQ(ierr); |
| } |
| /* Adjust values of D */ |
| ierr = VecGetArray(r,&pr);CHKERRQ(ierr); |
| ierr = VecGetArray(p,&pp);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(r,&pr);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(p,&pp);CHKERRQ(ierr); |
| ierr = VecGetArray(eps->D,&pD);CHKERRQ(ierr); |
| for (i=0;i<eps->nloc;i++) { |
| if (eps->balance == EPS_BALANCE_TWOSIDE) { |
| if (pp[i]!=0.0) pD[i] *= 1.0/PetscAbsScalar(pp[i]); |
| } |
| } |
| ierr = VecRestoreArray(r,&pr);CHKERRQ(ierr); |
| ierr = VecRestoreArray(p,&pp);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(r,&pr);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(p,&pp);CHKERRQ(ierr); |
| ierr = VecRestoreArray(eps->D,&pD);CHKERRQ(ierr); |
| } |
| PetscScalar alpha, Vec *X, PetscInt cX, PetscInt dX, const PetscScalar *M, |
| PetscInt ldM, PetscInt rM, PetscInt cM) |
| { |
| PetscErrorCode ierr; |
| PetscScalar *px, *py; |
| PetscInt rX, rY, ldX, ldY, i, rcX; |
| PetscErrorCode ierr; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt rX, rY, ldX, ldY, i, rcX; |
| PetscFunctionBegin; |
| if (rX != rY) { |
| SETERRQ(((PetscObject)*Y)->comm,1, "The multivectors do not have the same dimension"); |
| } |
| ierr = VecGetArray(X[0], &px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(X[0], &px);CHKERRQ(ierr); |
| ierr = VecGetArray(Y[0], &py);CHKERRQ(ierr); |
| /* Update the strides */ |
| ierr = SlepcDenseMatProd(py, ldY, beta, alpha, px, ldX, rX, rcX, |
| PETSC_FALSE, M, ldM, rM, cM, PETSC_FALSE); CHKERRQ(ierr); |
| ierr = VecRestoreArray(X[0], &px);CHKERRQ(ierr); |
| ierr = PetscObjectStateDecrease((PetscObject)X[0]); CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(X[0], &px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(Y[0], &py);CHKERRQ(ierr); |
| for(i=1; i<cM; i++) { |
| ierr = PetscObjectStateIncrease((PetscObject)Y[dY*i]); CHKERRQ(ierr); |
| const PetscScalar *M, PetscInt ldM, PetscInt rM, PetscInt cM, |
| PetscScalar *work, PetscInt lwork) |
| { |
| PetscErrorCode ierr; |
| PetscScalar **px, *Y, *Z; |
| PetscInt rX, i, j, rY, rY0, ldY; |
| PetscErrorCode ierr; |
| PetscScalar **px, *Y, *Z; |
| PetscInt rX, i, j, rY, rY0, ldY; |
| PetscFunctionBegin; |
| Vec *V, PetscInt sV, PetscInt eV, |
| PetscScalar *workS0, PetscScalar *workS1) |
| { |
| PetscErrorCode ierr; |
| PetscInt ldU, ldV, i, j, k; |
| PetscScalar *pu, *pv, *W, *Wr; |
| PetscErrorCode ierr; |
| PetscInt ldU, ldV, i, j, k; |
| const PetscScalar *pu, *pv; |
| PetscScalar *W, *Wr; |
| PetscFunctionBegin; |
| if (ldU != ldV) { |
| SETERRQ(((PetscObject)*U)->comm,1, "Matrix dimensions do not match"); |
| } |
| ierr = VecGetArray(U[0], &pu);CHKERRQ(ierr); |
| ierr = VecGetArray(V[0], &pv);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(U[0], &pu);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(V[0], &pv);CHKERRQ(ierr); |
| if (workS0) |
| W = workS0; |
| ierr = PetscFree(W); CHKERRQ(ierr); |
| } |
| ierr = VecRestoreArray(U[0], &pu); CHKERRQ(ierr); |
| ierr = PetscObjectStateDecrease((PetscObject)U[0]); CHKERRQ(ierr); |
| ierr = VecRestoreArray(V[0], &pv); CHKERRQ(ierr); |
| ierr = PetscObjectStateDecrease((PetscObject)V[0]); CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(U[0], &pu); CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(V[0], &pv); CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| Vec *V, PetscInt sV, PetscInt eV, DvdReduction *r, |
| DvdMult_copy_func *sr) |
| { |
| PetscErrorCode ierr; |
| PetscInt ldU, ldV; |
| PetscScalar *pu, *pv, *W; |
| PetscErrorCode ierr; |
| PetscInt ldU, ldV; |
| const PetscScalar *pu, *pv; |
| PetscScalar *W; |
| PetscFunctionBegin; |
| if (ldU != ldV) { |
| SETERRQ(((PetscObject)*U)->comm,1, "Matrix dimensions do not match"); |
| } |
| ierr = VecGetArray(U[0], &pu);CHKERRQ(ierr); |
| ierr = VecGetArray(V[0], &pv);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(U[0], &pu);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(V[0], &pv);CHKERRQ(ierr); |
| ierr = PetscLogEventBegin(SLEPC_VecsMult,0,0,0,0);CHKERRQ(ierr); |
| ierr = PetscLogEventEnd(SLEPC_VecsMult,0,0,0,0);CHKERRQ(ierr); |
| ierr = VecRestoreArray(U[0], &pu); CHKERRQ(ierr); |
| ierr = PetscObjectStateDecrease((PetscObject)U[0]); CHKERRQ(ierr); |
| ierr = VecRestoreArray(V[0], &pv); CHKERRQ(ierr); |
| ierr = PetscObjectStateDecrease((PetscObject)V[0]); CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(U[0], &pu); CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(V[0], &pv); CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "ShellMatMult_Cyclic" |
| static PetscErrorCode ShellMatMult_Cyclic(Mat B,Vec x,Vec y) |
| { |
| PetscErrorCode ierr; |
| SVD svd; |
| SVD_CYCLIC *cyclic; |
| PetscScalar *px,*py; |
| PetscInt m; |
| PetscErrorCode ierr; |
| SVD svd; |
| SVD_CYCLIC *cyclic; |
| const PetscScalar *px; |
| PetscScalar *py; |
| PetscInt m; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(B,(void**)&svd);CHKERRQ(ierr); |
| cyclic = (SVD_CYCLIC *)svd->data; |
| ierr = SVDMatGetLocalSize(svd,&m,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(cyclic->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(cyclic->x2,px+m);CHKERRQ(ierr); |
| ierr = VecResetArray(cyclic->x2);CHKERRQ(ierr); |
| ierr = VecResetArray(cyclic->y1);CHKERRQ(ierr); |
| ierr = VecResetArray(cyclic->y2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "SVDSetUp_Cyclic" |
| PetscErrorCode SVDSetUp_Cyclic(SVD svd) |
| { |
| PetscErrorCode ierr; |
| SVD_CYCLIC *cyclic = (SVD_CYCLIC *)svd->data; |
| PetscInt M,N,m,n,i,nloc,isl; |
| PetscScalar *pU,*isa,*va; |
| PetscBool trackall; |
| Vec v; |
| Mat Zm,Zn; |
| PetscErrorCode ierr; |
| SVD_CYCLIC *cyclic = (SVD_CYCLIC *)svd->data; |
| PetscInt M,N,m,n,i,nloc,isl; |
| const PetscScalar *isa; |
| PetscScalar *pU,*va; |
| PetscBool trackall; |
| Vec v; |
| Mat Zm,Zn; |
| PetscFunctionBegin; |
| ierr = MatDestroy(&cyclic->mat);CHKERRQ(ierr); |
| for (i=0; i<-svd->nini; i++) { |
| ierr = MatGetVecs(cyclic->mat,&v,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(v,&va);CHKERRQ(ierr); |
| ierr = VecGetArray(svd->IS[i],&isa);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(svd->IS[i],&isa);CHKERRQ(ierr); |
| ierr = VecGetSize(svd->IS[i],&isl);CHKERRQ(ierr); |
| if (isl == m) { |
| ierr = PetscMemcpy(va,isa,sizeof(PetscScalar)*m);CHKERRQ(ierr); |
| SETERRQ(((PetscObject)svd)->comm,PETSC_ERR_SUP,"Size of the initial subspace vectors should match to some dimension of A"); |
| } |
| ierr = VecRestoreArray(v,&va);CHKERRQ(ierr); |
| ierr = VecRestoreArray(svd->IS[i],&isa);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(svd->IS[i],&isa);CHKERRQ(ierr); |
| ierr = VecDestroy(&svd->IS[i]);CHKERRQ(ierr); |
| svd->IS[i] = v; |
| } |
| #define __FUNCT__ "SVDSolve_Cyclic" |
| PetscErrorCode SVDSolve_Cyclic(SVD svd) |
| { |
| PetscErrorCode ierr; |
| SVD_CYCLIC *cyclic = (SVD_CYCLIC *)svd->data; |
| PetscInt i,j,M,N,m,n; |
| PetscScalar sigma,*px; |
| Vec x,x1,x2; |
| PetscErrorCode ierr; |
| SVD_CYCLIC *cyclic = (SVD_CYCLIC *)svd->data; |
| PetscInt i,j,M,N,m,n; |
| PetscScalar sigma; |
| const PetscScalar *px; |
| Vec x,x1,x2; |
| PetscFunctionBegin; |
| ierr = EPSSolve(cyclic->eps);CHKERRQ(ierr); |
| ierr = EPSGetEigenpair(cyclic->eps,i,&sigma,PETSC_NULL,x,PETSC_NULL);CHKERRQ(ierr); |
| if (PetscRealPart(sigma) > 0.0) { |
| svd->sigma[j] = PetscRealPart(sigma); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(x2,px+m);CHKERRQ(ierr); |
| ierr = VecCopy(x1,svd->U[j]);CHKERRQ(ierr); |
| ierr = VecScale(svd->V[j],1.0/sqrt(2.0));CHKERRQ(ierr); |
| ierr = VecResetArray(x1);CHKERRQ(ierr); |
| ierr = VecResetArray(x2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| j++; |
| } |
| } |
| */ |
| PetscErrorCode MatLaplacian2D_Mult(Mat A,Vec x,Vec y) |
| { |
| void *ctx; |
| int nx,lo,j,one=1; |
| PetscScalar *px,*py,dmone=-1.0; |
| PetscErrorCode ierr; |
| void *ctx; |
| int nx,lo,j,one=1; |
| const PetscScalar *px; |
| PetscScalar *py,dmone=-1.0; |
| PetscErrorCode ierr; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(A,&ctx);CHKERRQ(ierr); |
| nx = *(int*)ctx; |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| tv(nx,&px[0],&py[0]); |
| tv(nx,&px[lo],&py[lo]); |
| BLASaxpy_(&nx,&dmone,&px[lo-nx],&one,&py[lo],&one); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |
| } |
| #define __FUNCT__ "MatBrussel_Mult" |
| PetscErrorCode MatBrussel_Mult(Mat A,Vec x,Vec y) |
| { |
| PetscInt n; |
| PetscScalar *px,*py; |
| CTX_BRUSSEL *ctx; |
| PetscErrorCode ierr; |
| PetscInt n; |
| const PetscScalar *px; |
| PetscScalar *py; |
| CTX_BRUSSEL *ctx; |
| PetscErrorCode ierr; |
| PetscFunctionBegin; |
| ierr = MatShellGetContext(A,(void**)&ctx);CHKERRQ(ierr); |
| ierr = MatGetLocalSize(ctx->T,&n,PETSC_NULL);CHKERRQ(ierr); |
| ierr = VecGetArray(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x1,px);CHKERRQ(ierr); |
| ierr = VecPlaceArray(ctx->x2,px+n);CHKERRQ(ierr); |
| ierr = VecAXPY(ctx->y2,-ctx->beta,ctx->x1);CHKERRQ(ierr); |
| ierr = VecAXPY(ctx->y2,-ctx->alpha * ctx->alpha + ctx->sigma,ctx->x2);CHKERRQ(ierr); |
| ierr = VecRestoreArray(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x1);CHKERRQ(ierr); |
| ierr = VecResetArray(ctx->x2);CHKERRQ(ierr); |
| @*/ |
| PetscErrorCode SlepcVecMAXPBY(Vec y,PetscScalar beta,PetscScalar alpha,PetscInt nv,PetscScalar a[],Vec x[]) |
| { |
| PetscErrorCode ierr; |
| PetscBLASInt n,m,one=1; |
| PetscScalar *py,*px; |
| PetscErrorCode ierr; |
| PetscBLASInt n,m,one=1; |
| PetscScalar *py; |
| const PetscScalar *px; |
| PetscFunctionBegin; |
| PetscValidHeaderSpecific(y,VEC_CLASSID,1); |
| ierr = PetscLogEventBegin(SLEPC_VecMAXPBY,*x,y,0,0);CHKERRQ(ierr); |
| ierr = VecGetArray(y,&py);CHKERRQ(ierr); |
| ierr = VecGetArray(*x,&px);CHKERRQ(ierr); |
| ierr = VecGetArrayRead(*x,&px);CHKERRQ(ierr); |
| n = PetscBLASIntCast(nv); |
| m = PetscBLASIntCast((y)->map->n); |
| BLASgemv_("N",&m,&n,&alpha,px,&m,a,&one,&beta,py,&one); |
| ierr = VecRestoreArray(y,&py);CHKERRQ(ierr); |
| ierr = VecRestoreArray(*x,&px);CHKERRQ(ierr); |
| ierr = VecRestoreArrayRead(*x,&px);CHKERRQ(ierr); |
| ierr = PetscLogFlops(nv*2*(y)->map->n);CHKERRQ(ierr); |
| ierr = PetscLogEventEnd(SLEPC_VecMAXPBY,*x,y,0,0);CHKERRQ(ierr); |
| PetscFunctionReturn(0); |