Seidel

/*
 * =====================================================================================
 *
 *       Filename:  exp_seidel.c
 *
 *    Description:  Solve example in page 74 - GT GTS LETRONGVINH
 *
 *        Version:  1.0
 *        Created:  06/09/2011 10:02:20 PM
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  FamiHug (fh), famihug@gmail.com
 *        Company:  familug.blogspot.com
 *
 * =====================================================================================
 */

//Todo List
//nhap ma tran
//viet ve dang alphaX + beta
//tach alpha thanh alpha 1 + alpha2
//thuc hien lap seidel

#include <stdio.h>
#include <string.h>
#include <stdlib.h>


//Nhap ma tran

void NhapMaTran(float *matrix)
{
    int size = 3;
    int i,j;

    printf("Nhap vao ma tran: \n");
    for(i = 0; i < 3; i++)
    {
        for(j = 0; j < 4; j++)
        {
            scanf("%f",matrix + i * (size + 1)  + j);
        }
    }
    printf("Nhap xong!\n");
    return;
}

void InMaTran(float *matrix)
{
    int size = 3;
    int i,j;

    printf("Ma Tran: \n");
    for(i = 0; i < size; i++)
    {
        for(j = 0; j < size + 1; j++)
        {
            printf("%3.3f ",*(matrix + i * (size + 1) + j));
        }
        printf("\n");
    }
    return;
}

//Viet lai theo dang alphaX
float *TraVeAlpha(float *matrix)
{
    int i,j;
    float *alpha;
    int size = 3;
    alpha = NULL;

    alpha = matrix;
    for(i = 0; i < size; i++)
    {
        for(j = 0; j < size; j++)
        {
            if(i != j) *(alpha + i * (size + 1) + j) = -(*(alpha + i * (size + 1) + j))/(*(alpha + i * (size + 1) + i));
        }
    }
    for(i = 0; i < 3; i++)
    {
        (*(alpha + i * (size + 1) + size)) /= (*(alpha + i * (size + 1) + i));
        *(alpha + i * (size + 1) + i) = 0;
    }
    return alpha;
}

//tach alpha thanh alpha1 + alpha 2
void TachAlpha(float *alpha, float *al1, float *al2)
{
    int i,j;
    int size = 3;

    //Copy alpha -> al1, al2
    memcpy(al1, alpha, sizeof(float) * size * (size + 1));
    memcpy(al2, alpha, sizeof(float) * size * (size + 1));

    //xoa nua tren duong cheo cua alpha de dc al1
    for(i = 0; i < size; i++)
    {
        for(j = 0; j < size; j++)
        {
            if(j >= i) *(al1 + i * (size + 1) + j) = 0;
            else *(al2 + i * (size + 1) + j) = 0;
        }
    }


    //Xoa nua duoi duong cheo cua matran alpha de dc al2


    return;
}


//Kiem tra tinh lap don

float absolute(float x)
{
    if(x >= 0) return x;
    else return -x;
}

int isLapDon(float **alpha)
{
    int i,j;
    int size = 3;
    float max = 0;
    float tempSum = 0;
    for(i = 0; i < 3; i++)
    {
        tempSum = 0;
        for(j = 0; j < 3; j++)
        {
            tempSum += absolute(*(*alpha + i * (size + 1) + j));
        }
        if(max < tempSum) max = tempSum;
    }
    printf("Max: %f\n",max);
    if(max < 1) return 1;
    else return 0;
}


//Chon xap xi dau, so lan lap k
//Thuc hien lap va dua ra ket qua

int main(int argc, char **argv)
{
    int size = 3;
    float *matrix;
    float *alpha, *al1, *al2;

    matrix = (float *)malloc(sizeof(float) * size * (size + 1));
    al1 = (float *)malloc(sizeof(float) * size * (size + 1));
    al2 = (float *)malloc(sizeof(float) * size * (size + 1));

    NhapMaTran(matrix);
    InMaTran(matrix);
    alpha = TraVeAlpha(matrix);
    InMaTran(alpha);
    printf("Co tinh lap don: %d\n",isLapDon(&alpha));

    TachAlpha(alpha,al1,al2);

    InMaTran(al1);
    InMaTran(al2);

    //Nhap vao xap xi dau
    int k = 4; //k la so lan lap
    int i,g;

    float R[3][10];
    R[0][0] = 0, R[1][0] = 0, R[2][0] = 0;

    for(g = 1; g < k + 1; g++)
    {
        for(i = 0; i < 3; i++)
        {
            R[i][g] = R[0][g - 1] * (*(al1 + i * (size + 1) + 0)) + R[1][g - 1] * (*(al1 + i * (size + 1) + 1)) + R[2][g - 1]  * (*(al1 + i * (size + 1) + 2)) + (*(al2 + i * (size + 1) + 0)) + R[1][g - 1] * (*(al2 + i * (size + 1) + 1)) + R[2][g - 1]  * (*(al2 + i * (size + 1) + 2)) + *(alpha + i * (size + 1) + size);
        }
    }

    printf("Ket Qua: \n");

    for(i = 0; i < 3; i++)
    {
        for(g = 0; g < k + 1; g++)
        {
            printf("%f ",R[i][g]);
        }
        printf("\n");
    }

    return 0;
}