Hello PETSc team,
firtst of all sorry if this is not the right email where I can address my
question to.
I was looking for an IDR(s) implementation on PETSc since it would be
needed for the research I was performing.
I did not find it (I don't know if it is on-working or not contemplated) so
I tried to implement it by myself (since I have Matlab code available)
Code works (solution is exact) but performance are not good as expected
(not so bad, but lightly higher than BicGSTAB)
I would like to share with you and, if someone has a little bit of time to
spend on, I would be grateful if he could help me
I share my preliminary code, sorry for bad "informatic" implementation, I'm
sure is much improvable from this point of view.
It works only in serial mode for the moment
Thank you all for the time
Regards
Simone Cammarasana
___________________________________
Saluti / Regards
*Simone Cammarasana*
Tel. +39 340 0029498
void idrSolve(Mat *A, Vec *b) {
PetscErrorCode ierr;
PetscInt one = 1.0, zero = 0.0;
PetscViewer viewer;
int i,j,k;
PetscInt sizeMat;
PetscScalar alfa = 1, beta = 0;
PetscInt s = 4;
PetscReal tol = 1e-12;
PetscReal angle = 0.7;
PetscReal om = 1;
PetscReal Alpha,Beta;
PetscInt iter = 0;
PetscInt maxIteration = 10000;
Mat P,G,U,M;
PetscScalar *aP, *aG, *aU, *aM, *ar, *af, *aGapp, *aUapp, *ax, *at; /* scalar from array */
PetscScalar *ac, *av;
Vec x0; /* initial solution */
Vec menox;
Vec r,f;
Vec D; /* Diagonal ones */
Vec Gapp, Uapp;
Vec t;
Vec v;
Vec c;
Vec x;
PetscReal normb;
PetscReal normr;
PetscReal tolb;
/*Read sizeMat */
ierr = VecGetSize(b,&sizeMat);
/*Create Matrix and Vectors */
ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,sizeMat,s,NULL,&P); ierr = MatSetFromOptions(P);
ierr = VecCreate(PETSC_COMM_WORLD, &x0); ierr = VecSetSizes(x0,PETSC_DECIDE,sizeMat); ierr = VecSetFromOptions(x0); ierr = VecSet(x0,zero);
ierr = VecCreate(PETSC_COMM_WORLD, &r); ierr = VecSetSizes(r,PETSC_DECIDE,sizeMat); ierr = VecSetFromOptions(r);
ierr = VecCreate(PETSC_COMM_WORLD, &menox); ierr = VecSetSizes(menox,PETSC_DECIDE,sizeMat); ierr = VecSetFromOptions(menox);
ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,sizeMat,s,NULL,&G); ierr = MatSetFromOptions(G);
ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,sizeMat,s,NULL,&U); ierr = MatSetFromOptions(U);
ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,s,s,NULL,&M); ierr = MatSetFromOptions(M);
ierr = VecCreate(PETSC_COMM_WORLD, &D); ierr = VecSetSizes(D,PETSC_DECIDE,s); ierr = VecSetFromOptions(D); ierr = VecSet(D,one);
ierr = VecCreate(PETSC_COMM_WORLD, &f); ierr = VecSetSizes(f,PETSC_DECIDE,s); ierr = VecSetFromOptions(f); ierr = VecSet(f,zero);
ierr = VecCreate(PETSC_COMM_WORLD, &Gapp); ierr = VecSetSizes(Gapp,PETSC_DECIDE,sizeMat); ierr = VecSetFromOptions(Gapp);
ierr = VecCreate(PETSC_COMM_WORLD, &Uapp); ierr = VecSetSizes(Uapp,PETSC_DECIDE,sizeMat); ierr = VecSetFromOptions(Uapp);
ierr = VecCreate(PETSC_COMM_WORLD, &t); ierr = VecSetSizes(t,PETSC_DECIDE,sizeMat); ierr = VecSetFromOptions(t);
ierr = VecCreate(PETSC_COMM_WORLD, &v); ierr = VecSetSizes(v,PETSC_DECIDE,sizeMat); ierr = VecSetFromOptions(v);
ierr = VecCreate(PETSC_COMM_WORLD, &c); ierr = VecSetSizes(c,PETSC_DECIDE,s); ierr = VecSetFromOptions(c);
ierr = VecCreate(PETSC_COMM_WORLD, &x); ierr = VecSetSizes(x,PETSC_DECIDE,sizeMat); ierr = VecSetFromOptions(x);
/*Create Random Object for P */
ierr = PetscMalloc1(sizeMat*s,&aP);
for(i = 0; i < sizeMat*s;i++) aP[i] = rand()/(double)RAND_MAX;
/*P Orthnorm to be performed here */
/* initialize x to x0 */
x = x0;
/* norm and relative tolerance */
ierr = VecNorm(b,NORM_2,&normb);
tolb = tol*normb;
/* Residual Calculation */
menox = x; /* surely to be optimize */
ierr = VecScale(menox, -1);
ierr = MatMultAdd(A,menox,b,r);
ierr = VecNorm(r,NORM_2,&normr);
/* Create M */
ierr = MatZeroEntries(M);
ierr = MatDiagonalSet(M,D,INSERT_VALUES);
/* Extract Arrays */
ierr = VecGetArray(r,&ar);
ierr = VecGetArray(f,&af);
ierr = VecGetArray(x,&ax);
ierr = MatDenseGetArray(M,&aM);
ierr = MatDenseGetArray(U,&aU);
ierr = MatDenseGetArray(G,&aG);
ierr = VecGetArray(Uapp,&aUapp);
ierr = VecGetArray(t,&at);
ierr = VecGetArray(v,&av);
ierr = VecGetArray(c,&ac);
while( (iter < maxIteration) && (normr > tolb) ) {
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
1, s, sizeMat, alfa, ar,1, aP, sizeMat, beta, af, 1); /*(m,n,k,alfa,A,lda,B,ldb,beta,C,ldc) */
for(k = 0; k < s; k++) {
/* Initialize ac */
for(i = k; i< s; i++) ac[i] = 0;
/* Solve M\f */
ac[k] = af[k]/aM[k*s+k];
for(i = k+1; i < s; i++) {
for(j = k; j < i; j++) {
ac[i] = ac[i]+(aM[j*s+i]*ac[j])/aM[i*s+i];
}
ac[i] = af[i]/aM[i*s+i]-ac[i];
}
/* Calculate v */
cblas_dcopy(sizeMat,ar,1,av,1);
cblas_dgemv(CblasColMajor, CblasNoTrans,
sizeMat, s-k, -1, aG+sizeMat*k,sizeMat, ac+k,1, 1, av, 1); /*layout, trans, m(A),n(A),alfa,A,lda,x,incx,beta,y,incy */
/* Compute new U */
cblas_dgemv(CblasColMajor, CblasNoTrans,
sizeMat, s-k, +1, aU+sizeMat*k,sizeMat, ac+k,1, om, av, 1); /*layout, trans, m(A),n(A),alfa,A,lda,x,incx,beta,y,incy --> y = alfa*Ax + beta*y */
cblas_dcopy(sizeMat,av,1,aU+sizeMat*k,1);
/*Compute new G */
cblas_dcopy(sizeMat,aU+k*sizeMat,1,aUapp,1);
ierr = VecRestoreArray(Uapp,&aUapp);
ierr = MatMult(A,Uapp,Gapp);
ierr = VecGetArray(Gapp,&aGapp);
ierr = VecGetArray(Uapp,&aUapp);
cblas_dcopy(sizeMat,aGapp,1,aG+sizeMat*k,1);
/*BiOrthogonalise the new basis vectors */
for(i = 0; i < k-1; i++) {
Alpha = cblas_ddot(sizeMat,aP+sizeMat*i,1,aG+sizeMat*k,1)/aM[i*s+i]; /*n,x,incx,y,incy */
cblas_daxpy(sizeMat,-Alpha,aG+i*sizeMat,1,aG+k*sizeMat,1);/* n,alfa,x,incx,y,incy --> y = ax+y */
cblas_daxpy(sizeMat,-Alpha,aU+i*sizeMat,1,aU+k*sizeMat,1);/* n,alfa,x,incx,y,incy --> y = ax+y */
}
/* New column of M */
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
1, s-k, sizeMat, alfa, aG+k*sizeMat,1, aP+k*sizeMat, sizeMat, 0, aM+k*s+k, 1); /*(m,n,k,alfa,A,lda,B,ldb,beta,C,ldc)-->C = alfaAB+betaC */
/* Make r ortoghonal */
Beta = af[k]/aM[k*s+k];
cblas_daxpy(sizeMat,-Beta,aG+k*sizeMat,1,ar,1);/* n,alfa,x,incx,y,incy --> y = ax+y */
cblas_daxpy(sizeMat,Beta,aU+k*sizeMat,1,ax,1);/* n,alfa,x,incx,y,incy --> y = ax+y */
normr = calcNorm(ar,sizeMat);
/* New f = P'*r */
if(k < s-1) cblas_daxpy(s-k-1,-Beta,aM+k*s+k+1,1,af+k+1,1); /*n,alfa,x,incx,y,incy --> y = ax+y*/
}
/*Calculate t = A*v */
cblas_dcopy(sizeMat,ar,1,av,1);
ierr = VecRestoreArray(v,&av);
ierr = MatMult(A,v,t);
ierr = VecGetArray(v,&av);
/*Calculate Omega */
om = omega(at,ar,angle,sizeMat);
/* Calculate r and x */
cblas_daxpy(sizeMat,-om,at,1,ar,1);/*n,alfa,x,incx,y,incy --> y = ax+y */
cblas_daxpy(sizeMat,+om,av,1,ax,1);/*n,alfa,x,incx,y,incy --> y = ax+y */
normr = calcNorm(ar,sizeMat);
iter++;
}
/* RESTORE */
ierr = VecRestoreArray(x,&ax);
/* STDOUT PRINT */
printf("iterations = %i\n",iter);
VecView(x,PETSC_VIEWER_STDOUT_WORLD);
return;
}
PetscReal omega(PetscScalar *at, PetscScalar *ar, PetscReal angle, PetscInt sizeMat) {
PetscReal ns,nt,ts;
PetscReal rho,om;
ns = calcNorm(ar,sizeMat);
nt = calcNorm(at,sizeMat);
ts = cblas_ddot(sizeMat,at,1,ar,1);
rho = fabs(ts/(nt*ns));
om = ts/(nt*nt);
if(rho < angle) om = om*angle/rho;
return om;
}
PetscReal calcNorm(PetscScalar *array, PetscInt size) {
int i;
PetscReal somma = 0;
for(i = 0; i < size; i++) somma = somma + array[i]*array[i];
return sqrt(somma);
}
