/*==============================================================================
  ==============================================================================
  CRIHAN     10/99
  Guy Moebs   Guy.Moebs@crihan.fr

  Utilisation de la bibliotheque scientifique BLAS (www.netlib.org) en C

  Calcul du vecteur R = ( U*(tU) + LAMBDA*V*(tV) )  *  ( U + V )

  FUNCTION(S) APPELEE(S) :
      DDOT    produit scalaire de 2 vecteurs                   en double precision 

  SUBROUTINE(S) APPELEE(S) :
      DGER    construction de la matrice A <- a * U * (tV) + A en double precision
      DAXPY   construction du vecteur Y <- Y + a * X           en double precision
      DGEMV   calcul du vecteur Y <- a * A * x + b * Y         en double precision

 ===============================================================================
 ===============================================================================

      DECLARATIONS

 ===============================================================================
 ===============================================================================*/

     /* INCLUDE(S) NECESSAIRE(S) */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <scsl_cblas.h>
#include <ccdim.h>
#include <ccscal.h>

     /* DECLARATION DES VARIABLES GLOBALES */

void main()
{
     /* DECLARATION DES VARIABLES LOCALES */

      double *a;
      double *u;
      double *v;

      double       lambda;
      double       alpha;
      double       beta;
      double       scal;

      int    incu;
      int    incv;
      int    i;
      int    j;

/*            MatrixTranspose  mtype = NoTranspose;  */
      const enum CBLAS_TRANSPOSE mtype = CblasNoTrans;
      const enum CBLAS_ORDER order = CblasColMajor;

      int    n  = NDIM;

     /* DECLARATION DES FUNCTION(S) EXTERIEURE(S) */

/*==============================================================================
  ==============================================================================

      DEBUT DU PROGRAMME

  ==============================================================================
  ==============================================================================*/

     /* Initialisation des vecteurs U et V */

      u = (double *) malloc ( n * sizeof(double) );
      v = (double *) malloc ( n * sizeof(double) );

      scal = (0.5e+1 * un) / 0.12e+2;
      for (i=0; i < n; i++)
      {
	 u [i] = scal;
         v [i] = 1.2e+0 + i;
      }

     /* Initialisation de la matrice A */

      a = (double *) malloc ( n*n * sizeof(double) );

      for ( i=0; i < n; i++ )
      {
	 for ( j=0; j < n; j++ )
	 {
	    a [n*i+j] = zero;
         }
      }

     /* Calcul de LAMBDA = U'*V */

      incu = 1;
      incv = 1;

      lambda = cblas_ddot ( n, u, incu, v, incv );


	/* Calcul de la matrice A = un * U * (tU) + A */

      alpha = un;

      cblas_dger ( order, n, n, alpha, u, incu, u, incu, a, n );


	/* Calcul de la matrice A = LAMBDA * V * (tV) + A */

      cblas_dger ( order, n, n, lambda, v, incv, v, incv, a, n );


	/* Calcul du vecteur V = U + V */

      alpha = un;

      cblas_daxpy ( n, alpha, u, incu, v, incv );


	/* Calcul du U = ALPHA * A * V + BETA * U */

      alpha = un;
      beta  = zero;

      cblas_dgemv (order, mtype, n, n, alpha, a, n, v, incv, beta, u, incu);


	/* Affichage des resultats */

      printf ("Resultats : %g, %g, %g, %g\n", u [0], u [1], u [n-2], u [n-1]);

      free (u);
      free (v);
      free (a);

/*==============================================================================
  ==============================================================================

     FIN DU PROGRAMME

  ==============================================================================
  ==============================================================================*/
      return;
}