CFD Online Discussion Forums

CFD Online Discussion Forums (
-   Main CFD Forum (
-   -   Problems with biCGSTAB (

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?

        r = addVectors(b, multiplyR(A, x0), n, -1);
        copyVector(r, rtilde, n);
        int i;
                rho_1 = dotProduct(rtilde, r, n);
                        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);
                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);
                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!!

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.