/* fichero maxvalp.c */
/* cálculo del máximo valor propio por itaración directa */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void main(void)
{
        int i, j,orden;
        double z, z2, x, x2;
        double **a, *vect, *vect2;
        double **ma_alloc(int, int);
        double vabs(double);
        void normalizar(double *, int);
        void matvec(double **, double*, double*, int);
        printf ("\nIntroduzca el orden de la matriz: ");
        scanf ("%d", &orden);
        a = ma_alloc(orden, orden);
        vect = malloc(orden*sizeof(double));
        vect2 = malloc(orden*sizeof(double));
        printf("\nIntroduzca los datos de la matriz:\n");
        for (i=0; i<orden; i++)
                for(j=0;j<orden;j++) {
                        printf("\nElemento (%d,%d): ", (i+1), (j+1));
                        scanf("%lf", &a[i][j]);
                }
        for (i=0; i<orden; i++)   /* aproximación inicial */
                vect2[i]=1.0;
        do {
                for (j=0; j<orden; j++)
                        vect[j] = vect2[j];
                normalizar (vect, orden);
                matvec(a, vect, vect2, orden);
                for (j=0; j<orden; j++)
                        vect[j] = vect2[j];
                for (j=0; j<orden; j++)
                        z += vect[j]*vect[j];
                x = sqrt (z);
                for (j=0; j<orden; j++)
                        z += vect2[j]*vect2[j];
                x2 = sqrt (z);
        } while (vabs((x2-x)/x2)<0.0000001);
        printf ("El vector buscado es el:\n(");
        for (i=0; i<orden; i++)
                printf ("%lf ", vect[i]);
        printf(")");
        printf("\nEl valor propio es: %lf", x);
}
double **ma_alloc(int filas, int columnas){
        int i;
        double **a;
        a = malloc(filas*sizeof(double *));
        for (i=0; i<filas; i++)
                a[i] = malloc(columnas*sizeof(double));
        return a;
}
void normalizar(double *v, int n){
        double x, z=0.0;
        int i;
        for (i=0; i<n; i++)
                z += v[i]*v[i];
        x = sqrt (z);
        for (i=0; i<n; i++)
                v[i] = v[i]/x;
}
void matvec(double **matriz, double *vector, double *vector2, int orden)
{
        int i,j;
        double sum;
        for (i=0; i<orden; i++) {
                sum = 0;
                for(j=0; j<orden; j++)
                        sum += matriz[i][j]*vector[j];
                vector2[i] = sum;
        }
}
double vabs(double a){
        if (a<0.0)
                return (-a);
        else
                return a;
}