CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   Problems with biCGSTAB (http://www.cfd-online.com/Forums/main/88274-problems-bicgstab.html)

 Petrov May 12, 2011 11:03

Problems with biCGSTAB

I wrote a biconjugate gradient method in C for my FEM solver, but it does not work. (I can tell because when I use a different method, the results are correct). Can you see if there is an obvious problem with the code below?
Code:

```        r = addVectors(b, multiplyR(A, x0), n, -1);         copyVector(r, rtilde, n);         int i;         for(i=0;i<n;i++)         {                 rho_1 = dotProduct(rtilde, r, n);                 if(i==0)                 {                         copyVector(r, p, n);                 } else {                         beta = (rho_1/rho_2) * (alpha/omega);                         double * tempA = addVectors(p, v, n, -omega);                         p = addVectors(r, tempA, n, beta);                         free(tempA);                 }                 v = multiplyR(A, p);                 alpha = rho_1 / (dotProduct(rtilde, v, n));                 s = addVectors(r, v, n, (-1)*alpha);                 t = multiplyR(A, s);                 omega = (dotProduct(t, s, n))/dotProduct(t, t, n);                 double* temp = addVectors(x, p, n, alpha);                 x = addVectors(temp, s, n, omega);                 free(temp);                 r = addVectors(s, t, n, -omega);                 rho_2 = rho_1;         }```
The code you see is the "meat" of the solver. Before comes the allocation, and after comes the free-ing of memory. The input of the function is A, the matrix; b the vector, and x0 the initial guess.

multiplyR(matrix, vector) is a function that multiplies a matrix on the left with a vector on the right.

addVectors(vector A, vector B, int n, double k) returns A+kB. (The n is length of each vector).

As you can see, I'm not using a pre-conditioner.

Any help would be much appreciated!!

:::UPDATE:::
I have recently tried using a gauss-siedel preconditioner. The results did not improve at all, but they did change. Does that mean that I should try more sophisticated preconditioners, that there something wrong with the algorithm above, or that I should move on to gmres or some other method?

(I think I'm just going to try gaussian elimination until I feel like tackling this problem again. I've spent far too long on it!!!)

 All times are GMT -4. The time now is 09:08.