/*------------------------------------------------------------------------------*
 * File Name: OCN_f06.h															*
 * Creation: TCZ 7/25/2001														*
 * Purpose: Origin C Header for NAG functions									*
 * Copyright (c); OriginLab Corp.	2001										*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *------------------------------------------------------------------------------*/

#ifndef _O_NAG_f06_H
#define _O_NAG_f06_H

//#importdll "ONAG" // NAG DLL prepared by OriginLab
#pragma dll(ONAG)
#include <NAG\nag_types.h>
/* 
////////////////////////////////////////////////////////////////////////////////
This headfile contains the basic linear algebra functions in NAG_Chapter f06 which 
perform elementary algebraic operations involving vectors and matricces. 

The enumeration types given by
  typedef enum 
  { 
      NoTranspose,  //Operate with the matrix
      Transpose,    //Operate with the transpose of the matrix
      ConjugateTranspose   //Operate with the conjugate transpose of the matrix
  } MatrixTranspose;
  
  typedef enum 
  { 
      UpperTriangle, //Upper triangle
      LowerTriangle  // Lower triangle
  } MatrixTriangle;
  
  typedef enum 
  { 
      UnitTriangular, // Unit triangular, the diagonal elements are not referenced
      NotUnitTriangular // Non-unit triangular
  } MatrixUnitTriangular;
  
  typedef enum 
  { 
      LeftSide, // Multiply general matrix by symmetric, Hermitian or triangular matrix on the left
      RightSide // Multiply general matrix by symmetric, Hermitian or triangular matrix on the right
  } OperationSide;
*/
/////////////////////////////////////////////////////////////////////////////////

/** f06pac
		computes a matrix-vector product for general matrix.
		
Example:
	This example computes the matrix-vector product for a general matrix
        1  2  3  4
        5  6  7  8
        9 10 11 12
       13 14 15 16
    
    and a vector
       17
       18
       19
       20.

void test_f06pac()      //test for function dgemv(f06pac)
{
    int nRows=4;
    int nCols=4;

    matrix<double> mat_input={{1,2,3,4},{5,6,7,8},{9,10,11,12},{13,14,15,16}};

    vector<double> vec_input={17,18,19,20};

    vector<double> vec_output(nRows);
    vec_output=0;

    dgemv(NoTranspose,nRows,nCols,1,mat_input,nCols,vec_input,1,1,vec_output,1);

    printf("Compute a matrix-vector product for general real matrix M*V:\n");

    printf("M=\n");                     //output mat_input
    for(int i=0;i<nRows;i++)
    {
       printf("       ");
       for(int j=0;j<nCols;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_input
    printf("V=\n");
    for(i=0;i<nRows;i++)
    {
      printf("       ");
      printf("%f\n",vec_input[i]);
    }

    printf("\n");                       //output vec_output
    printf("Result vector=\n");
    for(i=0;i<nRows;i++)
    {
      printf("       ");
      printf("%f\n",vec_output[i]);
    }
}

	The output is following:
	
Compute a matrix-vector product for general real matrix M*V:
M=
       1.000000    2.000000    3.000000    4.000000    
       5.000000    6.000000    7.000000    8.000000    
       9.000000    10.000000    11.000000    12.000000    
       13.000000    14.000000    15.000000    16.000000    

V=
       17.000000
       18.000000
       19.000000
       20.000000

Result vector=
       190.000000
       486.000000
       782.000000
       1078.000000


Parameters:
	trans=NoTranspose, return (alpha*Ax+beta*y) in array y[];
	trans=Transpose, return (alpha*transpose(A)*x+beta*y) in array y[];
	trans=ConjugateTranspose, return (alpha*ConjugateTranspose(A)*x+beta*y) in array y[].
Return: NULL
*/

void  dgemv(
MatrixTranspose trans,
int m,  //the number of rows of A
int n,  // the number of columns of A
double alpha, 
const double a[], //a general m by n matrix
int tda, //the last demension of A
const double x[], //vector x
int incx,  //the spacing which the elements of the vector x occur
double beta,
double y[], //Input: vector y; Output: alpha*trans(A)*x + beta*y
int incy   //the spacing which the elements of the vector y occur
);

/** f06sac
       Compute a matrix-vector product for complex general matrix
Example:
	This example computes the matrix-vector product for a general complex matrix
        1+2i,8+3i,2+3i,3+4i
        2+1i,3+1i,5+2i,4+4i
        1+3i,2+1i,2+3i,5+4i
        1+4i,2+3i,6+1i,2+4i
    
    and a vector
       2+2i
       1+1i
       3+2i
       2+4i.

void test_f06sac()    //test for function zgemv(f06sac)
{
   int nRows=4;
   int nCols=4;

   matrix<complex>  mat_input={{1+2i,8+3i,2+3i,3+4i},
                               {2+1i,3+1i,5+2i,4+4i},
                               {1+3i,2+1i,2+3i,5+4i},
                               {1+4i,2+3i,6+1i,2+4i}};

   vector<complex>  vec_input={2+2i,1+1i,3+2i,2+4i};

   vector<complex>  vec_output(nRows);
   vec_output=0;

   zgemv(NoTranspose,nRows,nCols,1,mat_input,nCols,vec_input,1,1,vec_output,1);

   printf("Compute a matrix-vector product for general complex matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRows;i++)
   {
       printf("       ");
       for(int j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}

	The output is following:
	
Compute a matrix-vector product for general complex matrix M*V:
M=
       (1.000000,2.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (2.000000,1.000000i)    (3.000000,1.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (1.000000,3.000000i)    (2.000000,1.000000i)    (2.000000,3.000000i)    (5.000000,4.000000i)    
       (1.000000,4.000000i)    (2.000000,3.000000i)    (6.000000,1.000000i)    (2.000000,4.000000i)    

V=
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)

Result vector=
       (-7.000000,50.000000i)
       (7.000000,50.000000i)
       (-9.000000,52.000000i)
       (-3.000000,46.000000i)

Parameters:
	trans=NoTranspose, return (alpha*Ax+beta*y) in array y[];
	trans=Transpose, return (alpha*transpose(A)*x+beta*y) in array y[];
	trans=ConjugateTranspose, return (alpha*ConjugateTranspose(A)*x+beta*y) in array y[].
Return: NULL
*/
void  zgemv(
MatrixTranspose trans, 
int m,  //the number of rows of A
int n,  //the number of columns of A
complex alpha, 
const complex a[], //the m by n complex matrix A 
int tda, // the last dimension of A
const complex x[], // the vector x
int incx,  //the spacing which the elements of the vector x occur
complex beta, 
complex y[], //Input: the vecor y; Output: alpha*trans(A)*x+beta*y
int incy   // the spacing which the elements of the vector y occur
);



/** f06pbc
       Compute a matrix-vector product for real general band matrix
Example:
	This example computes the matrix-vector product for a real matrix
	   1.000000    2.000000    3.000000    0.000000    
       4.000000    5.000000    6.000000    7.000000    
       0.000000    8.000000    9.000000    10.000000  
    which in band storage and a real vector
       11 
       12
       13
       14
       
void test_f06pbc()              //test for function dgbmv(f06pbc)
{
   int nRows=3;
   int nCols=4;

   matrix<double> mat_input={{1,2,3,0},
                             {4,5,6,7},
                             {0,8,9,10}};
                             
   matrix<double> mat_input_band={{0,1,2,3},
                                  {4,5,6,7},
                                  {8,9,10,0}};
   int kupper=2;
   int klower=1;

   vector<double> vec_input={11,12,13,14};

   vector<double> vec_output(nRows);
   vec_output=0;

   dgbmv(NoTranspose,nRows,nCols,klower,kupper,1,mat_input_band,nCols,vec_input,1,1,vec_output,1);

   printf("Compute a matrix-vector product for real band matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRows;i++)
   {
       printf("       ");
       for(int j=0;j<nCols;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nCols;i++)
   {
      printf("       ");
      printf("%f\n",vec_input[i]);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRows;i++)
   {
      printf("       ");
      printf("%f\n",vec_output[i]);
   }
}

  The output is following:
Compute a matrix-vector product for real band matrix M*V:
M=
       1.000000    2.000000    3.000000    0.000000    
       4.000000    5.000000    6.000000    7.000000    
       0.000000    8.000000    9.000000    10.000000    

V=
       11.000000
       12.000000
       13.000000
       14.000000

Result vector=
       74.000000
       280.000000
       353.000000
  	

Parameters:
	trans=NoTranspose, return (alpha*Ax+beta*y) in array y[];
	trans=Transpose, return (alpha*transpose(A)*x+beta*y) in array y[];
	trans=ConjugateTranspose, return (alpha*ConjugateTranspose(A)*x+beta*y) in array y[].
Return: NULL
*/
void  dgbmv(
MatrixTranspose trans, 
int m, // the number of rows of A
int n, // the number of columns of A 
int kl, //specify the subdiagonals of the matrix A
int ku, //specify the superdiagonals of the matrix A
double alpha,
const double a[], //the real m by n matrix A which stored in band storage
int tda, // the last dimension of A
const double x[], // vector x
int incx,  // the spacing which the elements of the vector x occur
double beta, 
double y[], //Input: vector y; Output: alpha*trans(A)*x+beta*y
int incy   // the spacing which the elements of vector y occur
);


/** f06sbc
       Compute a matrix-vector product for complex general band matrix
Example:
	This example computes the matrix-vector product for a complex general matrix
	   (1.000000,1.200000i)    (2.000000,2.100000i)    (3.000000,1.000000i)    (0.000000,0.000000i)    
       (4.000000,2.600000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (7.000000,4.100000i)    
       (0.000000,0.000000i)    (8.000000,1.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
    which is stored in band storage and a complex vector
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)
       
void test_f06sbc()     //test for function zgbmv(f06sbc)
{
   int nRows=3;
   int nCols=4;

   matrix<complex> mat_input={{1+1.2i,2+2.1i,3+1.0i,0},
                             {4+2.6i,5+1.2i,6+2.6i,7+4.1i},
                             {0,8+1i,9,10+2i}};
   int kupper=2;
   int klower=1;
   matrix<complex> mat_input_band={{0,1+1.2i,2+2.1i,3+1.0i},
                                  {4+2.6i,5+1.2i,6+2.6i,7+4.1i},
                                  {8+1i,9,10+2i,0}};

   vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};

   vector<complex> vec_output(nRows);
   vec_output=0;

   zgbmv(NoTranspose,nRows,nCols,klower,kupper,1,mat_input_band,nCols,vec_input,1,1,vec_output,1);

   printf("Compute a matrix-vector product for complex band matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRows;i++)
   {
       printf("       ");
       for(int j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nCols;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}

   The output is as following:
Compute a matrix-vector product for complex band matrix M*V:
M=
       (1.000000,1.200000i)    (2.000000,2.100000i)    (3.000000,1.000000i)    (0.000000,0.000000i)    
       (4.000000,2.600000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (7.000000,4.100000i)    
       (0.000000,0.000000i)    (8.000000,1.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    

V=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (69.600000,59.400000i)
       (262.220000,166.800000i)
       (349.400000,76.000000i)
   	
Parameters:
	trans=NoTranspose, return (alpha*Ax+beta*y) in array y[];
	trans=Transpose, return (alpha*transpose(A)*x+beta*y) in array y[];
	trans=ConjugateTranspose, return (alpha*ConjugateTranspose(A)*x+beta*y) in array y[].
Return: NULL
*/

void  zgbmv(
MatrixTranspose trans, 
int m, //the number of rows of A
int n,  // the number of columns of A
int kl, //the subdiagonals of matrix A
int ku, // the superdiagonals of matrix A
complex alpha,
const complex a[], //the m by n general complex matrix A
int tda, //the last dimension of A
const complex x[], //complex vector x
int incx, //the spacing which the elements of the vector x occur
complex beta, 
complex y[], //Input:complex vector y; Output: alpha*trans(A)*x+beta*y
int incy  // the spacing which the elements of the vector y occur
);



/** f06pcc
       Compute a matrix-vector product for real symmetric matrix
Example:
	This example computes the matrix-vector product for a real symmetric matrix
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       2.000000    7.000000    9.000000    10.000000    
       4.000000    8.000000    10.000000    11.000000    

    and a real vector
       12.000000
       13.000000
       14.000000
       15.000000

       
void test_f06pcc()         //test for function dsymv(f06pcc)
{
    int nColRow=4;

    matrix<double> mat_input={{1,2,3,4},
                              {2,6,7,8},
                              {2,7,9,10},
                              {4,8,10,11}};

    vector<double> vec_input={12,13,14,15};

    vector<double> vec_output(nColRow);
    vec_output=0;

    dsymv(UpperTriangle,nColRow,1,mat_input,nColRow,vec_input,1,1,vec_output,1);

    printf("Compute a matrix-vector product for real symmetric matrix M*V:\n");

    printf("M=\n");                     //output mat_input
    for(int i=0;i<nColRow;i++)
    {
       printf("       ");
       for(int j=0;j<nColRow;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_input
    printf("V=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_input[i]);
    }

    printf("\n");                       //output vec_output
    printf("Result vector=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_output[i]);
    }
}

   The output is as following:
Compute a matrix-vector product for real symmetric matrix M*V:
M=
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       2.000000    7.000000    9.000000    10.000000    
       4.000000    8.000000    10.000000    11.000000    

V=
       12.000000
       13.000000
       14.000000
       15.000000

Result vector=
       140.000000
       320.000000
       403.000000
       457.000000

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.
Return: NULL        

*/
void  dsymv(
MatrixTriangle uplo, 
int n,  //the order of symmetric matrix A 
double alpha,
const double a[], //the real symmetric matrix of order n
int tda,  //the last dimension of A
const double x[], // vector x
int incx,  //the spacing which the elements of vector x occur
double beta, 
double y[],  //Input: vector y; Output: alpha*Ax+beta*y
int incy
);


/** f06scc
       Compute a matrix-vector product for complex Hermitian matrix
Example:
	This example computes the matrix-vector product for a complex Hermitian matrix
       (1.000000,0.000000i)    (2.000000,2.100000i)    (3.000000,1.000000i)    (4.000000,1.300000i)    
       (2.000000,-2.100000i)   (5.000000,0.000000i)    (6.000000,2.600000i)    (7.000000,4.100000i)    
       (3.000000,-1.000000i)   (6.000000,-2.600000i)   (9.000000,0.000000i)    (10.000000,2.000000i)    
       (4.000000,-1.300000i)   (7.000000,-4.100000i)   (10.000000,-2.000000i)  (11.000000,0.000000i)    

    and a complex vector
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)


       
void test_f06scc()      //test for function zhemv(f06scc)
{
   int nRowCol=4;

   matrix<complex> mat_input={{1,2+2.1i,3+1.0i,4+1.3i},
                             {2-2.1i,5,6+2.6i,7+4.1i},
                             {3-1.0i,6-2.6i,9,10+2i},
                             {4-1.3i,7-4.1i,10-2i,11}};

   vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};

   vector<complex> vec_output(nRowCol);
   vec_output=0;

   zhemv(UpperTriangle,nRowCol,1,mat_input,nRowCol,vec_input,1,1,vec_output,1);

   printf("Compute a matrix-vector product for complex Hermitian matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}

   The output is as following:
Compute a matrix-vector product for complex Hermitian matrix M*V:
M=
       (1.000000,0.000000i)    (2.000000,2.100000i)    (3.000000,1.000000i)    (4.000000,1.300000i)    
       (2.000000,-2.100000i)    (5.000000,0.000000i)    (6.000000,2.600000i)    (7.000000,4.100000i)    
       (3.000000,-1.000000i)    (6.000000,-2.600000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (4.000000,-1.300000i)    (7.000000,-4.100000i)    (10.000000,-2.000000i)    (11.000000,0.000000i)    

V=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (125.660000,71.600000i)
       (249.620000,96.700000i)
       (360.400000,27.800000i)
       (418.600000,-41.700000i)

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.
Return: NULL        
*/
void  zhemv(
MatrixTriangle uplo, 
int n, // the order of complex Hermitian matrix A
complex alpha,
const complex a[], // order n complex Hermitian matrix A
int tda, // the last dimension of A
const complex x[], //complex vector
int incx, //the spacing which the elements of the vector x occur
complex beta, 
complex y[], //Input: complex vector y; Output: alpha*Ax+beta*y 
int incy  // the spacing which the elements of the vector y occuur
);

/** f06pdc
       Compute a matrix-vector product for symmetric band matrix
Example:
	This example computes the matrix-vector product for a symmetric band matrix
       1.000000    2.000000    3.000000    0.000000    
       2.000000    5.000000    6.000000    7.000000    
       3.000000    6.000000    9.000000    10.000000    
       0.000000    7.000000    10.000000   11.000000    

    and a real vector
       11.000000
       12.000000
       13.000000
       14.000000
       

void test_f06pdc()              //test for function dsbmv(f06pdc)
{
   int nRowCol=4;

   matrix<double> mat_input={{1,2,3,0},
                             {2,5,6,7},
                             {3,6,9,10},
                             {0,7,10,11}};
   int k_nzero=2;
   matrix<double> mat_input_band={{1,2,3},
                                  {5,6,7},
                                  {9,10,0},
                                  {11,0,0}};

   vector<double> vec_input={11,12,13,14};

   vector<double> vec_output(nRowCol);
   vec_output=0;

   dsbmv(UpperTriangle,nRowCol,k_nzero,1,mat_input_band,k_nzero+1,vec_input,1,1,vec_output,1);

   printf("Compute a matrix-vector product for real symmetric band matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("%f\n",vec_input[i]);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("%f\n",vec_output[i]);
   }
}

   The output is as following:
Compute a matrix-vector product for real symmetric band matrix M*V:
M=
       1.000000    2.000000    3.000000    0.000000    
       2.000000    5.000000    6.000000    7.000000    
       3.000000    6.000000    9.000000    10.000000    
       0.000000    7.000000    10.000000    11.000000    

V=
       11.000000
       12.000000
       13.000000
       14.000000

Result vector=
       74.000000
       258.000000
       362.000000
       368.000000

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.
Return: NULL
*/

void  dsbmv(
MatrixTriangle uplo, 
int n, // the order of symmetric band matrix A
int k, // the subdiagonals (superdiagonals) of A
double alpha, 
const double a[], //the symmetric band matrix A
int tda, //the last dimension of A
const double x[], //vector x
int incx, //the spacing which the elements of vector x occur
double beta,
double y[], //Input: vector y; Output: alpha*Ax+beta*y
int incy //the spacing which the elements of vector y occur
);

/** f06sdc
       Compute a matrix-vector product for complex Hermitian band matrix
Example:
	This example computes the matrix-vector product for a complex Hermitian band matrix
       (1.000000,0.000000i)    (2.000000,2.100000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    
       (2.000000,-2.100000i)    (5.000000,0.000000i)    (6.000000,2.600000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (6.000000,-2.600000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (10.000000,-2.000000i)    (11.000000,0.000000i)    


    and a complex vector
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)
       

void test_f06sdc()      //test for function zhbmv(f06sdc)
{
   int nRowCol=4;

   matrix<complex> mat_input={{1,2+2.1i,0,0},
                             {2-2.1i,5,6+2.6i,0},
                             {0,6-2.6i,9,10+2i},
                             {0,0,10-2i,11}};

   int k_nzero=1;
   matrix<complex> mat_input_band={{1,2+2.1i},
                                   {5,6+2.6i},
                                   {9,10+2i},
                                   {11,0}}; 
                             
   vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};

   vector<complex> vec_output(nRowCol);
   vec_output=0;

   zhbmv(UpperTriangle,nRowCol,k_nzero,1,mat_input_band,k_nzero+1,vec_input,1,1,vec_output,1);

   printf("Compute a matrix-vector product for complex Hermitian band matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}

   The output is as following:
Compute a matrix-vector product for complex Hermitian band matrix M*V:
M=
       (1.000000,0.000000i)    (2.000000,2.100000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    
       (2.000000,-2.100000i)    (5.000000,0.000000i)    (6.000000,2.600000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (6.000000,-2.600000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (10.000000,-2.000000i)    (11.000000,0.000000i)    

V=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (35.000000,27.200000i)
       (159.000000,26.700000i)
       (325.400000,32.800000i)
       (288.000000,13.800000i)


Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.
Return: NULL
*/

void  zhbmv(
MatrixTriangle uplo, 
int n, //the order of Hermitian band matrix A
int k, //the subdiagonals (superdiagonals) of A
complex alpha, 
const complex a[], //the complex Hermitian band matrix A
int tda, //the last dimension of A
const complex x[], //vector x
int incx,  //the spacing which the elements of vector x occur
complex beta,
complex y[], //Input: vector y; Output: alpha*Ax+beta*y
int incy  // the spacing which the elements of vector y occur
);


/** f06pec
       Compute a matrix-vector product for real symmetric packed matrix
Example:
	This example computes the matrix-vector product for a real symmetric packed matrix
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       2.000000    7.000000    9.000000    10.000000    
       4.000000    8.000000    10.000000    11.000000  

    and a real vector
       12.000000
       13.000000
       14.000000
       15.000000

void test_f06pec()         //test for function dspmv(f06pec)
{
    int nColRow=4;

    matrix<double> mat_input_pack={1,2,3,4,
                                     6,7,8,
                                       9,10,
                                         11};

    matrix<double> mat_input={{1,2,3,4},
                              {2,6,7,8},
                              {2,7,9,10},
                              {4,8,10,11}};

    vector<double> vec_input={12,13,14,15};

    vector<double> vec_output(nColRow);
    vec_output=0;

    dspmv(UpperTriangle,nColRow,1,mat_input_pack,vec_input,1,1,vec_output,1);

    printf("Compute a matrix-vector product for real symmetric packed matrix M*V:\n");

    printf("M=\n");                     //output mat_input
    for(int i=0;i<nColRow;i++)
    {
       printf("       ");
       for(int j=0;j<nColRow;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_input
    printf("V=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_input[i]);
    }

    printf("\n");                       //output vec_output
    printf("Result vector=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_output[i]);
    }
}

   The output is as following:
Compute a matrix-vector product for real symmetric packed matrix M*V:
M=
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       2.000000    7.000000    9.000000    10.000000    
       4.000000    8.000000    10.000000    11.000000    

V=
       12.000000
       13.000000
       14.000000
       15.000000

Result vector=
       140.000000
       320.000000
       403.000000
       457.000000

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.
Return: NULL
*/
void  dspmv(
MatrixTriangle uplo, 
int n, // the order of the real symmetric packed matrix A 
double alpha, 
const double ap[], //the real symmetric matrix A in packed storage 
const double x[], //vetor x
int incx,  // the spacing which the elements of vector x occur
double beta, 
double y[], //Input: vector y; Output: alpha*Ax+beta*y
int incy  // the spacing which the elements of vector y occur
);

/** f06sec
       Compute a matrix-vector product for complex Hermitian packed matrix
Example:
	This example computes the matrix-vector product for a complex Hermitian packed matrix
       (1.000000,0.000000i)     (2.000000,2.100000i)    (3.000000,1.000000i)    (4.000000,1.300000i)    
       (2.000000,-2.100000i)    (5.000000,0.000000i)    (6.000000,2.600000i)    (7.000000,4.100000i)    
       (3.000000,-1.000000i)    (6.000000,-2.600000i)   (9.000000,0.000000i)    (10.000000,2.000000i)    
       (4.000000,-1.300000i)    (7.000000,-4.100000i)   (10.000000,-2.000000i)  (11.000000,0.000000i)    

    and a complex vector
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)


void test_f06sec()      //test for function zhpmv(f06sec)
{
   int nRowCol=4;

   matrix<complex> mat_input={{1,2+2.1i,3+1.0i,4+1.3i},
                             {2-2.1i,5,6+2.6i,7+4.1i},
                             {3-1.0i,6-2.6i,9,10+2i},
                             {4-1.3i,7-4.1i,10-2i,11}};
       
   matrix<complex> mat_input_pack={1,2+2.1i,3+1.0i,4+1.3i,
                                          5,6+2.6i,7+4.1i,
                                                       9,10+2i,
                                                         11};

   vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};

   vector<complex> vec_output(nRowCol);
   vec_output=0;

   zhpmv(UpperTriangle,nRowCol,1,mat_input_pack,vec_input,1,1,vec_output,1);

   printf("Compute a matrix-vector product for complex packed Hermitian matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}

   The output is as following:
Compute a matrix-vector product for complex packed Hermitian matrix M*V:
M=
       (1.000000,0.000000i)    (2.000000,2.100000i)    (3.000000,1.000000i)    (4.000000,1.300000i)    
       (2.000000,-2.100000i)    (5.000000,0.000000i)    (6.000000,2.600000i)    (7.000000,4.100000i)    
       (3.000000,-1.000000i)    (6.000000,-2.600000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (4.000000,-1.300000i)    (7.000000,-4.100000i)    (10.000000,-2.000000i)    (11.000000,0.000000i)    

V=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (125.660000,71.600000i)
       (249.620000,96.700000i)
       (360.400000,27.800000i)
       (418.600000,-41.700000i)


Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.
Return: NULL
*/
void  zhpmv(
MatrixTriangle uplo, 
int n, //the order of the complex Hermitian matrix A
complex alpha, 
const complex ap[], // order n complex Hermitian matrix A in packed storage
const complex x[], //vector x
int incx, //the spacing which the elements of vector x occur
complex beta, 
complex y[], //Input: vector y; Output: alpha*Ax+beta*y
int incy //the spacing which the elements of vector y occur
);


/** f06pfc
       Compute a matrix-vector product for real triangular matrix
Example:
	This example computes the matrix-vector product for a real triangular matrix
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000      

    and a real vector
       12.000000
       13.000000
       14.000000
       15.000000

void test_f06pfc()         //test for function dtrmv(f06pfc)
{
    int nColRow=4;

    matrix<double> mat_input={{1,2,3,4},
                              {0,6,7,8},
                              {0,0,9,10},
                              {0,0,0,11}};

    vector<double> vec_input={12,13,14,15};

    vector<double> vec_output;
    vec_output=vec_input;
    
    dtrmv(UpperTriangle,NoTranspose,NotUnitTriangular,nColRow,
            mat_input,nColRow,vec_output,1);

    printf("Compute a matrix-vector product for real triangular matrix M*V:\n");

    printf("M=\n");                     //output mat_input
    for(int i=0;i<nColRow;i++)
    {
       printf("       ");
       for(int j=0;j<nColRow;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_input
    printf("V=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_input[i]);
    }

    printf("\n");                       //output vec_output
    printf("Result vector=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_output[i]);
    }
}

   The output is as following:
Compute a matrix-vector product for real triangular matrix M*V:
M=
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000    

V=
       12.000000
       13.000000
       14.000000
       15.000000

Result vector=
       140.000000
       296.000000
       276.000000
       165.000000

Parameters:
	trans=NoTranspose, return A*x in array x[];
	trans=Transpose, return transpose(A)*x in array x[];
	trans=ConjugateTranspose, return ConjugateTranspose(A)*x in array x[].
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.
Return: NULL
*/
void  dtrmv(
MatrixTriangle uplo, 
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n,  //the order of the real triangular matrix A
const double a[], // order n real triangular matrix A
int tda, // the last dimension of A
double x[], //Input: the vector x; Output: the product of A and x
int incx  // the spacing which the elements of x occur
);

/** f06sfc
       Compute a matrix-vector product for complex triangular matrix
Example:
	This example computes the matrix-vector product for a complex triangular matrix
       (1.000000,1.200000i)    (2.000000,2.100000i)    (4.000000,0.500000i)    (1.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (2.000000,1.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    

    and a complex vector
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)


void test_f06sfc()      //test for function ztrmv(f06sfc)
{
   int nRowCol=4;

   matrix<complex> mat_input={{1+1.2i,2+2.1i,4+0.5i,1},
                              {0,5+1.2i,6+2.6i,2+1i},
                               {0,0,9,10+2i},
                               {0,0,0,11+4i}};

   int k_nzero=1;

   vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};

   vector<complex> vec_output(nRowCol);
   vec_output=vec_input;

   ztrmv(UpperTriangle,NoTranspose,NotUnitTriangular,nRowCol,
            mat_input,nRowCol,vec_output,1);

   printf("Compute a matrix-vector product for complex Hermitian band matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}

   The output is as following:
Compute a matrix-vector product for complex Hermitian band matrix M*V:
M=
       (1.000000,1.200000i)    (2.000000,2.100000i)    (4.000000,0.500000i)    (1.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (2.000000,1.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    

V=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (97.600000,56.700000i)
       (159.000000,77.800000i)
       (253.400000,64.000000i)
       (146.800000,75.800000i)


Parameters:
	trans=NoTranspose, return A*x in array x[];
	trans=Transpose, return transpose(A)*x in array x[];
	trans=ConjugateTranspose, return ConjugateTranspose(A)*x in array x[].
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/

void  ztrmv(
MatrixTriangle uplo, 
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n,  //the order of the complex triangular matrix A
const complex a[], //order n complex triangular matrix A
int tda, //the last dimension of A
complex x[], // Input: the complex vector x; Output: the product of A and x
int incx //the spacing which the elemnets of x occur
);

/** f06pgc
       Compute a matrix-vector product for real triangular band matrix
Example:
	This example computes the matrix-vector product for a real triangular band matrix
       1.000000    2.000000    3.000000    0.000000    
       0.000000    5.000000    6.000000    7.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000    
  

    and a real vector
       11.000000
       12.000000
       13.000000
       14.000000

void test_f06pgc()              //test for function dtbmv(f06pgc)
{
   int nRowCol=4;

   matrix<double> mat_input={{1,2,3,0},
                             {0,5,6,7},
                             {0,0,9,10},
                             {0,0,0,11}};
   int k_nzero=2;
   matrix<double> mat_input_band={{1,2,3},
                                  {5,6,7},
                                  {9,10,0},
                                  {11,0,0}};
                             
   vector<double> vec_input={11,12,13,14};

   vector<double> vec_output(nRowCol);
   vec_output=vec_input;

   dtbmv(UpperTriangle,NoTranspose,NotUnitTriangular,nRowCol,k_nzero,
            mat_input_band,k_nzero+1,vec_output,1);

   printf("Compute a matrix-vector product for real tringular band matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("%f\n",vec_input[i]);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("%f\n",vec_output[i]);
   }
}

   The output is as following:
Compute a matrix-vector product for real tringular band matrix M*V:
M=
       1.000000    2.000000    3.000000    0.000000    
       0.000000    5.000000    6.000000    7.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000    

V=
       11.000000
       12.000000
       13.000000
       14.000000

Result vector=
       74.000000
       236.000000
       257.000000
       154.000000

Parameters:
	trans=NoTranspose, return A*x in array x[];
	trans=Transpose, return transpose(A)*x in array x[];
	trans=ConjugateTranspose, return ConjugateTranspose(A)*x in array x[].
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/

void  dtbmv(
MatrixTriangle uplo, 
MatrixTranspose trans,
MatrixUnitTriangular diag, 
int n, //the order of A
int k, //the subdiagonals(superdiagonals) of A
const double a[], //the real triangular band matrix A
int tda, //the last dimension of A
double x[], //Input: the vector x; Output: the product of A and x
int incx //the spacing which the elements of vector x occur
);


/** f06sgc
       Compute a matrix-vector product for complex triangular band matrix
Example:
	This example computes the matrix-vector product for a complex triangular band matrix
       (1.000000,1.200000i)    (2.000000,2.100000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    

    and a complex vector
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)


void test_f06sgc()      //test for function ztbmv(f06sgc)
{
   int nRowCol=4;

   matrix<complex> mat_input={{1+1.2i,2+2.1i,0,0},
                             {0,5+1.2i,6+2.6i,0},
                             {0,0,9,10+2i},
                             {0,0,0,11+4i}};

   int k_nzero=1;
   matrix<complex> mat_input_band={{1+1.2i,2+2.1i},
                             {5+1.2i,6+2.6i},
                             {9,10+2i},
                             {11+4i,0}};
                             
   vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};

   vector<complex> vec_output(nRowCol);
   vec_output=vec_input;

   ztbmv(UpperTriangle,NoTranspose,NotUnitTriangular,nRowCol,k_nzero,
             mat_input_band,k_nzero+1,vec_output,1);

   printf("Compute a matrix-vector product for complex triangular band matrix M*V:\n");

   printf("M=\n");                     //output mat_input
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}

   The output is as following:
Compute a matrix-vector product for complex triangular band matrix M*V:
M=
       (1.000000,1.200000i)    (2.000000,2.100000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    

V=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (32.600000,40.400000i)
       (132.800000,60.200000i)
       (253.400000,64.000000i)
       (146.800000,75.800000i)

Parameters:
	trans=NoTranspose, return A*x in array x[];
	trans=Transpose, return transpose(A)*x in array x[];
	trans=ConjugateTranspose, return ConjugateTranspose(A)*x in array x[].
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/

void  ztbmv(
MatrixTriangle uplo, 
MatrixTranspose trans,
MatrixUnitTriangular diag, 
int n, 
int k,
const complex a[], 
int tda, 
complex x[],
int incx
);

/** f06phc
       Compute a matrix-vector product for real triangular packed matrix
Example:
	This example computes the matrix-vector product for a real triangular packed matrix
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000 
       
    and a real vector
       12.000000
       13.000000
       14.000000
       15.000000


void test_f06phc()         //test for function dtpmv(f06phc)
{
    int nColRow=4;

    matrix<double> mat_input_pack={1,2,3,4,
                                     6,7,8,
                                       9,10,
                                         11};

    matrix<double> mat_input={{1,2,3,4},
                              {0,6,7,8},
                              {0,0,9,10},
                              {0,0,0,11}};

    vector<double> vec_input={12,13,14,15};

    vector<double> vec_output(nColRow);
    vec_output=vec_input;

    dtpmv(UpperTriangle,NoTranspose,NotUnitTriangular,nColRow,
             mat_input_pack,vec_output,1);

    printf("Compute a matrix-vector product for real tringular packed matrix M*V:\n");

    printf("M=\n");                     //output mat_input
    for(int i=0;i<nColRow;i++)
    {
       printf("       ");
       for(int j=0;j<nColRow;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_input
    printf("V=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_input[i]);
    }

    printf("\n");                       //output vec_output
    printf("Result vector=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_output[i]);
    }
}


   The output is as following:
Compute a matrix-vector product for real tringular packed matrix M*V:
M=
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000    

V=
       12.000000
       13.000000
       14.000000
       15.000000

Result vector=
       140.000000
       296.000000
       276.000000
       165.000000

Parameters:
	trans=NoTranspose, return A*x in array x[];
	trans=Transpose, return transpose(A)*x in array x[];
	trans=ConjugateTranspose, return ConjugateTranspose(A)*x in array x[].
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/
void  dtpmv(
MatrixTriangle uplo, 
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n, //the order of A
const double ap[], //the real triangular packed matrix A
double x[], //Input: the vector x; Output: the product of A and x
int incx //the spacing which the elements of x occur
);

/** f06shc
       Compute a matrix-vector product for complex triangular packed matrix
Example:
	This example computes the matrix-vector product for a complex triangular packed matrix
       (1.000000,1.200000i)    (2.000000,2.100000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    

    and a complex vector
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)



void test_f06shc()         //test for function ztpmv(f06shc)
{
    int nRowCol=4;

    matrix<complex> mat_input_pack={1+1.2i,2+2.1i,0,0,
                                      5+1.2i,6+2.6i,0,
                                              9,10+2i,
                                                11+4i};

    matrix<complex> mat_input={{1+1.2i,2+2.1i,0,0},
                               {0,5+1.2i,6+2.6i,0},
                               {0,0,9,10+2i},
                               {0,0,0,11+4i}};

    vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};

    vector<complex> vec_output(nRowCol);
    vec_output=vec_input;

    ztpmv(UpperTriangle,NoTranspose,NotUnitTriangular,nRowCol,
             mat_input_pack,vec_output,1);

    printf("Compute a matrix-vector product for complex tringular packed matrix M*V:\n");

    printf("M=\n");                     //output mat_input
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_input
   printf("V=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_output
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}


   The output is as following:
Compute a matrix-vector product for complex tringular packed matrix M*V:
M=
       (1.000000,1.200000i)    (2.000000,2.100000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    

V=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (32.600000,40.400000i)
       (132.800000,60.200000i)
       (253.400000,64.000000i)
       (146.800000,75.800000i)

Parameters:
	trans=NoTranspose, return A*x in array x[];
	trans=Transpose, return transpose(A)*x in array x[];
	trans=ConjugateTranspose, return ConjugateTranspose(A)*x in array x[].
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/
void  ztpmv(
MatrixTriangle uplo, 
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n, //the order of A
const complex ap[], //the complex triangular packed matrix A 
complex x[], //Input: the vector x; Output: the product of A and x
int incx //the spacing which the elements of x occur
);

/** f06pjc
       Solve a system of equations for real triangular coefficient matrix
Example:
	This example solve a system of equations Ax=b, where matrix A is
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000  

    and the vector b is
       12.000000
       13.000000
       14.000000
       15.000000


void test_f06pjc()         //test for function dtrsv(f06pjc)
{
    int nColRow=4;

    matrix<double> mat_input={{1,2,3,4},
                              {0,6,7,8},
                              {0,0,9,10},
                              {0,0,0,11}};

    vector<double> vec_input={12,13,14,15};

    vector<double> vec_output;
    vec_output=vec_input;
    
    dtrsv(UpperTriangle,NoTranspose,NotUnitTriangular,nColRow,
            mat_input,nColRow,vec_output,1);

    printf("Solve a system of equations for real triangular matrix M*X=B:\n");

    printf("M=\n");                     //output mat_inputM
    for(int i=0;i<nColRow;i++)
    {
       printf("       ");
       for(int j=0;j<nColRow;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_inputB
    printf("B=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_input[i]);
    }

    printf("\n");                       //output vec_outputX
    printf("Result vector=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_output[i]);
    }
}


   The output is as following:
Solve a system of equations for real triangular matrix M*X=B:
M=
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000    

B=
       12.000000
       13.000000
       14.000000
       15.000000

Result vector=
       5.821549
       0.301347
       0.040404
       1.363636

Parameters:
	trans=NoTranspose, solve the equations A*x=b;
	trans=Transpose, solve the equations Transpose(A)*x=b;
	trans=ConjugateTranspose, solve the equations ConjugateTranspose(A)*x=b.
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/
void  dtrsv(
MatrixTriangle uplo, 
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n, //the order of A
const double a[], //the order n real triangular coefficient matrix
int tda, //the last dimension of A
double x[],//Input: the vector b; Output: the solution of the equations, vector x
int incx // the spacing which the elements of vecter x occur
);

/** f06sjc
       Solve a system of equations for complex triangular coefficient matrix
Example:
	This example solve a system of equations Ax=b, where matrix A is
       (1.000000,1.200000i)    (2.000000,2.100000i)    (4.000000,0.500000i)    (1.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (2.000000,1.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    


    and the vector b is
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

void test_f06sjc()      //test for function ztrsv(f06sjc)
{
   int nRowCol=4;

   matrix<complex> mat_input={{1+1.2i,2+2.1i,4+0.5i,1},
                              {0,5+1.2i,6+2.6i,2+1i},
                               {0,0,9,10+2i},
                               {0,0,0,11+4i}};

   vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};

   vector<complex> vec_output(nRowCol);
   vec_output=vec_input;

   ztrsv(UpperTriangle,NoTranspose,NotUnitTriangular,nRowCol,
            mat_input,nRowCol,vec_output,1);

   printf("Solve a system of equations for complex triangular matrix M*X=B:\n");

   printf("M=\n");                     //output mat_inputM
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_inputB
   printf("B=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_outputX
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}


   The output is as following:
Solve a system of equations for complex triangular matrix M*X=B:
M=
       (1.000000,1.200000i)    (2.000000,2.100000i)    (4.000000,0.500000i)    (1.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (2.000000,1.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    

B=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (1.498443,-2.397366i)
       (1.702533,-0.884189i)
       (0.078345,0.254339i)
       (1.176642,-0.264234i)

Parameters:
	trans=NoTranspose, solve the equations A*x=b;
	trans=Transpose, solve the equations Transpose(A)*x=b;
	trans=ConjugateTranspose, solve the equations ConjugateTranspose(A)*x=b.
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/
void  ztrsv(
MatrixTriangle uplo, 
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n, // the order of A
const complex a[], //the order n complex triangular coefficient matrix A
int tda, // the last dimension of A
complex x[], //Input: the vctor b; Output: the solution of the equations, vector x
int incx //the spacing which the elements of vector x occur
);

/** f06pkc
       Solve a system of equations for real triangular band coefficient matrix
Example:
	This example solve a system of equations Ax=b, where matrix A is
       1.000000    2.000000    0.000000    0.000000    
       0.000000    1.000000    2.000000    0.000000    
       0.000000    0.000000    1.000000    2.000000    
       0.000000    0.000000    0.000000    1.000000   
       
    and the vector b is
       2.000000
       2.000000
       2.000000
       2.000000

void test_f06pkc()              //test for function dtbsv(f06pkc)
{
   int nRowCol=4;

   matrix<double> mat_input={{1,2,0,0},
                             {0,1,2,0},
                             {0,0,1,2},
                             {0,0,0,1}};
   int k_nzero=1;
   matrix<double> mat_input_band={{1,2},
                                 {1,2},
                                 {1,2},
                                 {1,0}};

   vector<double> vec_input={2,2,2,2};

   vector<double> vec_output(nRowCol);
   vec_output=vec_input;

   dtbsv(UpperTriangle,NoTranspose,NotUnitTriangular,nRowCol,k_nzero,
            mat_input_band,k_nzero+1,vec_output,1);
            

   printf("Solve a system of equations for real triangular band matrix M*X=B:\n");

   printf("M=\n");                     //output mat_inputM
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
   }

   printf("\n");                       //output vec_inputB
   printf("B=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("%f\n",vec_input[i]);
   }

   printf("\n");                       //output vec_outputX
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("%f\n",vec_output[i]);
   }
   
}


   The output is as following:
Solve a system of equations for real triangular band matrix M*X=B:
M=
       1.000000    2.000000    0.000000    0.000000    
       0.000000    1.000000    2.000000    0.000000    
       0.000000    0.000000    1.000000    2.000000    
       0.000000    0.000000    0.000000    1.000000    

B=
       2.000000
       2.000000
       2.000000
       2.000000

Result vector=
       -10.000000
       6.000000
       -2.000000
       2.000000

Parameters:
	trans=NoTranspose, solve the equations A*x=b;
	trans=Transpose, solve the equations Transpose(A)*x=b;
	trans=ConjugateTranspose, solve the equations ConjugateTranspose(A)*x=b.
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/
void  dtbsv(
MatrixTriangle uplo,
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n, //the oder of A
int k, //the subdiagonals(superdiagonals) of A
const double a[], //the real triangular band coefficient matrix A
int tda, //the last dimension of A
double x[], //Input: the vector b; Output: the solution of the equations, vector x
int incx //the spacing which the elements of vector x occur
);

/** f06skc
       Solve a system of equations for complex triangular band coefficient matrix
Example:
	This example solve a system of equations Ax=b, where matrix A is
       (1.000000,1.200000i)    (2.000000,2.100000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    
       
    and the vector b is
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

void test_f06skc()      //test for function ztbsv(f06skc)
{
   int nRowCol=4;

   matrix<complex> mat_input={{1+1.2i,2+2.1i,0,0},
                              {0,5+1.2i,6+2.6i,0},
                               {0,0,9,10+2i},
                               {0,0,0,11+4i}};

   vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};
   
   int k_nonzero=1;
   matrix<complex> mat_input_band={{1+1.2i,2+2.1i},
                                   {5+1.2i,6+2.6i},
                                   {9,10+2i},
                                   {11+4i,0}};

   vector<complex> vec_output(nRowCol);
   vec_output=vec_input;

   ztbsv(UpperTriangle,NoTranspose,NotUnitTriangular,nRowCol,k_nonzero,
            mat_input_band,k_nonzero+1,vec_output,1);

   printf("Solve a system of equations for complex triangular band matrix M*X=B:\n");

   printf("M=\n");                     //output mat_inputM
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_inputB
   printf("B=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_outputX
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}


   The output is as following:
Solve a system of equations for complex triangular band matrix M*X=B:
M=
       (1.000000,1.200000i)    (2.000000,2.100000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (0.000000,0.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    

B=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (1.474731,-2.685431i)
       (2.226943,-0.880413i)
       (0.078345,0.254339i)
       (1.176642,-0.264234i)

Parameters:
	trans=NoTranspose, solve the equations A*x=b;
	trans=Transpose, solve the equations Transpose(A)*x=b;
	trans=ConjugateTranspose, solve the equations ConjugateTranspose(A)*x=b.
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/
void  ztbsv(
MatrixTriangle uplo,
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n, // the order of A
int k, //the subdiagonals(superdiagonals) of A
const complex a[], //the complex triangular band coefficient matrix A
int tda, //the last dimension of A
complex x[], //Input: the vector b; Output: the solution of the equations, vector x
int incx //the spacing which the elements of vector x occur
);

/** f06plc
       Solve a system of equations for real triangular packed coefficient matrix
Example:
	This example solve a system of equations Ax=b, where matrix A is
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000   
       
    and the vector b is
       12.000000
       13.000000
       14.000000
       15.000000

void test_f06plc()         //test for function dtpsv(f06plc)
{
    int nColRow=4;

    matrix<double> mat_input={{1,2,3,4},
                              {0,6,7,8},
                              {0,0,9,10},
                              {0,0,0,11}};
                              
   matrix<double> mat_input_pack={1,2,3,4,
                                    6,7,8,
                                     9,10,
                                       11};
                             

    vector<double> vec_input={12,13,14,15};

    vector<double> vec_output;
    vec_output=vec_input;
    
    dtpsv(UpperTriangle,NoTranspose,NotUnitTriangular,nColRow,
            mat_input_pack,vec_output,1);

    printf("Solve a system of equations for real packed triangular matrix M*X=B:\n");

    printf("M=\n");                     //output mat_inputM
    for(int i=0;i<nColRow;i++)
    {
       printf("       ");
       for(int j=0;j<nColRow;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_inputB
    printf("B=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_input[i]);
    }

    printf("\n");                       //output vec_outputX
    printf("Result vector=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_output[i]);
    }
}


   The output is as following:
Solve a system of equations for real packed triangular matrix M*X=B:
M=
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    9.000000    10.000000    
       0.000000    0.000000    0.000000    11.000000    

B=
       12.000000
       13.000000
       14.000000
       15.000000

Result vector=
       5.821549
       0.301347
       0.040404
       1.363636

Parameters:
	trans=NoTranspose, solve the equations A*x=b;
	trans=Transpose, solve the equations Transpose(A)*x=b;
	trans=ConjugateTranspose, solve the equations ConjugateTranspose(A)*x=b.
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/

void  dtpsv(
MatrixTriangle uplo,
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n, //the order of A
const double ap[], //the real triangular packed coefficient matrix A
double x[], //Input: the vector b; Output: the solutions of the equations, vector x
int incx  //the spacing which the elements of vector x occur
);

/** f06slc
       Solve a system of equations for complex triangular packed coefficient matrix
Example:
	This example solve a system of equations Ax=b, where matrix A is
       (1.000000,1.200000i)    (2.000000,2.100000i)    (4.000000,0.500000i)    (1.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (2.000000,1.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    
       
    and the vector b is
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)


void test_f06slc()      //test for function ztpsv(f06slc)
{
   int nRowCol=4;

   matrix<complex> mat_input={{1+1.2i,2+2.1i,4+0.5i,1},
                              {0,5+1.2i,6+2.6i,2+1i},
                               {0,0,9,10+2i},
                               {0,0,0,11+4i}};

  matrix<complex> mat_input_pack={1+1.2i,2+2.1i,4+0.5i,1,
                                      5+1.2i,6+2.6i,2+1i,
                                                 9,10+2i,
                                                   11+4i};
                                                   
   vector<complex> vec_input={11+2i,12,13+2i,14+1.8i};

   vector<complex> vec_output(nRowCol);
   vec_output=vec_input;

   ztpsv(UpperTriangle,NoTranspose,NotUnitTriangular,nRowCol,
            mat_input_pack,vec_output,1);

   printf("Solve a system of equations for complex packed triangular matrix M*X=B:\n");

   printf("M=\n");                     //output mat_inputM
   for(int i=0;i<nRowCol;i++)
   {
       printf("       ");
       for(int j=0;j<nRowCol;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_inputB
   printf("B=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_input[i].m_re,vec_input[i].m_im);
   }

   printf("\n");                       //output vec_outputX
   printf("Result vector=\n");
   for(i=0;i<nRowCol;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_output[i].m_re,vec_output[i].m_im);
   }
}


   The output is as following:
Solve a system of equations for complex packed triangular matrix M*X=B:
M=
       (1.000000,1.200000i)    (2.000000,2.100000i)    (4.000000,0.500000i)    (1.000000,0.000000i)    
       (0.000000,0.000000i)    (5.000000,1.200000i)    (6.000000,2.600000i)    (2.000000,1.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (9.000000,0.000000i)    (10.000000,2.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (11.000000,4.000000i)    

B=
       (11.000000,2.000000i)
       (12.000000,0.000000i)
       (13.000000,2.000000i)
       (14.000000,1.800000i)

Result vector=
       (1.498443,-2.397366i)
       (1.702533,-0.884189i)
       (0.078345,0.254339i)
       (1.176642,-0.264234i)

Parameters:
	trans=NoTranspose, solve the equations A*x=b;
	trans=Transpose, solve the equations Transpose(A)*x=b;
	trans=ConjugateTranspose, solve the equations ConjugateTranspose(A)*x=b.
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
    diag specify whether or not the matrix is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/
void  ztpsv(
MatrixTriangle uplo,
MatrixTranspose trans, 
MatrixUnitTriangular diag, 
int n, // the order of A
const complex ap[], //the complex triangular packed coefficient matrix A
complex x[], //Input: the vector b; Output: the solution of the equations, vector x
int incx //the spacing which the elements of vector x occur
);

/** f06pmc
       Perform a rank-one update for real general matrix:A=alpha*x*transpose(y)+A
Example:
	This example perform a rank-one update for A, where matrix A is
       1.000000    2.000000    3.000000    4.000000    
       5.000000    6.000000    7.000000    8.000000    
       9.000000    10.000000    11.000000    12.000000    
       13.000000    14.000000    15.000000    16.000000    
       
    and the vector x is
       17.000000
       18.000000
       19.000000
       20.000000
     
    vector y is
       1.000000
       2.000000
       3.000000
       4.000000


void test_f06pmc()      //test for function dger(f06pmc)
{
    int nRows=4;
    int nCols=4;

    matrix<double> mat_input={{1,2,3,4},
                              {5,6,7,8},
                              {9,10,11,12},
                              {13,14,15,16}};

    vector<double> vec_x={17,18,19,20};
    vector<double> vec_y={1,2,3,4};

    matrix<double> mat_output;
    mat_output=mat_input;
    
    int i,j;

    dger(nRows,nCols,1,vec_x,1,vec_y,1,mat_output,nRows);

    printf("Perform a rank-one update for real general matrix:A=A+X*transpose(Y)\n");

    printf("A=\n");                     //output mat_inputA
    for(i=0;i<nRows;i++)
    {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_inputX
    printf("X=\n");
    for(i=0;i<nRows;i++)
    {
      printf("       ");
      printf("%f\n",vec_x[i]);
    }

    
    printf("\n");                       //output vec_inputY
    printf("Y=\n");
    for(i=0;i<nCols;i++)
    {
      printf("       ");
      printf("%f\n",vec_y[i]);
    }
    
    printf("Result matrix A=\n");       //output mat_outputA
    for(i=0;i<nRows;i++)
    {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("%f    ",mat_output[i][j]);
       printf("\n");
    }
}


   The output is as following:
Perform a rank-one update for real general matrix:A=A+X*transpose(Y)
A=
       1.000000    2.000000    3.000000    4.000000    
       5.000000    6.000000    7.000000    8.000000    
       9.000000    10.000000    11.000000    12.000000    
       13.000000    14.000000    15.000000    16.000000    

X=
       17.000000
       18.000000
       19.000000
       20.000000

Y=
       1.000000
       2.000000
       3.000000
       4.000000
Result matrix A=
       18.000000    36.000000    54.000000    72.000000    
       23.000000    42.000000    61.000000    80.000000    
       28.000000    48.000000    68.000000    88.000000    
       33.000000    54.000000    75.000000    96.000000    

Parameters:
Return: NULL
*/
void  dger(
int m, //the number of rows of A
int n, //the number of column of A
double alpha, 
const double x[], //the vector x 
int incx, //the spacing which the elements of x occur
const double y[], //the vector y
int incy, // the spacing which the elements of y occur
double a[], //Input: the real general matrix A; Output: the matrix A after rank-one update
int tda // the last dimension of A
);

/** f06smc
       Perform a rank-one update for complex general matrix:A=alpha*x*transpose(y)+A
Example:
	This example perform a rank-one update for A, where matrix A is
       (1.000000,2.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (2.000000,1.000000i)    (3.000000,1.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (1.000000,3.000000i)    (2.000000,1.000000i)    (2.000000,3.000000i)    (5.000000,4.000000i)    
       (1.000000,4.000000i)    (2.000000,3.000000i)    (6.000000,1.000000i)    (2.000000,4.000000i)    
   
    and the vector x is
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)

     
    vector y is
       (1.000000,2.000000i)
       (3.000000,0.400000i)
       (1.000000,0.200000i)
       (4.000000,2.100000i)

void test_f06smc()    //test for function zgeru(f06smc)
{
   int nRows=4;
   int nCols=4;

   matrix<complex>  mat_input={{1+2i,8+3i,2+3i,3+4i},
                               {2+1i,3+1i,5+2i,4+4i},
                               {1+3i,2+1i,2+3i,5+4i},
                               {1+4i,2+3i,6+1i,2+4i}};
                               

   vector<complex>  vec_x={2+2i,1+1i,3+2i,2+4i};
   vector<complex>  vec_y={1+2i,3+0.4i,1+0.2i,4+2.1i};

   matrix<complex>  mat_output();
   mat_output=mat_input;
   
   int i,j;

   zgeru(nRows,nCols,1,vec_x,1,vec_y,1,mat_output,nRows);

   printf("Perform a rank-one update for complex general matrix:A=A+X*transpose(Y)\n");

   printf("A=\n");                     //output mat_inputA
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_inputX
   printf("X=\n");
   for(i=0;i<nRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_x[i].m_re,vec_x[i].m_im);
   }

   printf("\n");                       //output vec_outputY
   printf("Y=\n");
   for(i=0;i<nCols;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_y[i].m_re,vec_y[i].m_im);
   }
   
   printf("Result matrix A=\n");                     //output mat_inputA
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_output[i][j].m_re,mat_output[i][j].m_im);
       printf("\n");
   }
}


   The output is as following:
Perform a rank-one update for complex general matrix:A=A+X*transpose(Y)
A=
       (1.000000,2.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (2.000000,1.000000i)    (3.000000,1.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (1.000000,3.000000i)    (2.000000,1.000000i)    (2.000000,3.000000i)    (5.000000,4.000000i)    
       (1.000000,4.000000i)    (2.000000,3.000000i)    (6.000000,1.000000i)    (2.000000,4.000000i)    

X=
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)

Y=
       (1.000000,2.000000i)
       (3.000000,0.400000i)
       (1.000000,0.200000i)
       (4.000000,2.100000i)
Result matrix A=
       (-1.000000,8.000000i)    (13.200000,9.800000i)    (3.600000,5.400000i)    (6.800000,16.200000i)    
       (1.000000,4.000000i)    (5.600000,4.400000i)    (5.800000,3.200000i)    (5.900000,10.100000i)    
       (0.000000,11.000000i)    (10.200000,8.200000i)    (4.600000,5.600000i)    (12.800000,18.300000i)    
       (-5.000000,12.000000i)    (6.400000,15.800000i)    (7.200000,5.400000i)    (1.600000,24.200000i)    


Parameters:
Return: NULL
*/
void  zgeru(
int m, //the number of rows of A
int n, //the number of columns of A
complex alpha, 
const complex x[], //the vecor x
int incx, //the spacing which the elements of x occur
const complex y[],//vector y
int incy, //the spacing which the elements of y occur
complex a[], //Input: the complex general matrix A; Output: matrix A after update
int tda //the last dimension of A
);


/** f06snc
       Perform a rank-one update for complex general matrix:A=alpha*x*ConjugateTranspose(y)+A
Example:
	This example perform a rank-one update for A, where matrix A is
       (1.000000,2.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (2.000000,1.000000i)    (3.000000,1.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (1.000000,3.000000i)    (2.000000,1.000000i)    (2.000000,3.000000i)    (5.000000,4.000000i)    
       (1.000000,4.000000i)    (2.000000,3.000000i)    (6.000000,1.000000i)    (2.000000,4.000000i)    
   
    and the vector x is
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)
     
    vector y is
       (1.000000,2.000000i)
       (3.000000,0.400000i)
       (1.000000,0.200000i)
       (4.000000,2.100000i)


void test_f06snc()    //test for function zgerc(f06snc)
{
   int nRows=4;
   int nCols=4;

   matrix<complex>  mat_input={{1+2i,8+3i,2+3i,3+4i},
                               {2+1i,3+1i,5+2i,4+4i},
                               {1+3i,2+1i,2+3i,5+4i},
                               {1+4i,2+3i,6+1i,2+4i}};
                               

   vector<complex>  vec_x={2+2i,1+1i,3+2i,2+4i};
   vector<complex>  vec_y={1+2i,3+0.4i,1+0.2i,4+2.1i};

   matrix<complex>  mat_output();
   mat_output=mat_input;
   
   int i,j;

   zgerc(nRows,nCols,1,vec_x,1,vec_y,1,mat_output,nRows);

   printf("Perform a rank-one update for complex general matrix:A=A+X*transpose(Y)\n");

   printf("A=\n");                     //output mat_inputA
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_inputX
   printf("X=\n");
   for(i=0;i<nRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_x[i].m_re,vec_x[i].m_im);
   }

   printf("\n");                       //output vec_outputY
   printf("Y=\n");
   for(i=0;i<nCols;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_y[i].m_re,vec_y[i].m_im);
   }
   
   printf("Result matrix A=\n");                     //output mat_inputA
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_output[i][j].m_re,mat_output[i][j].m_im);
       printf("\n");
   }
}


   The output is as following:
Perform a rank-one update for complex general matrix:A=A+X*transpose(Y)
A=
       (1.000000,2.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (2.000000,1.000000i)    (3.000000,1.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (1.000000,3.000000i)    (2.000000,1.000000i)    (2.000000,3.000000i)    (5.000000,4.000000i)    
       (1.000000,4.000000i)    (2.000000,3.000000i)    (6.000000,1.000000i)    (2.000000,4.000000i)    

X=
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)

Y=
       (1.000000,2.000000i)
       (3.000000,0.400000i)
       (1.000000,0.200000i)
       (4.000000,2.100000i)
Result matrix A=
       (7.000000,0.000000i)    (14.800000,8.200000i)    (4.400000,4.600000i)    (15.200000,7.800000i)    
       (5.000000,0.000000i)    (6.400000,3.600000i)    (6.200000,2.800000i)    (10.100000,5.900000i)    
       (8.000000,-1.000000i)    (11.800000,5.800000i)    (5.400000,4.400000i)    (21.200000,5.700000i)    
       (11.000000,4.000000i)    (9.600000,14.200000i)    (8.800000,4.600000i)    (18.400000,15.800000i)    


Parameters:
Return: NULL
*/
void  zgerc(
int m, //the number of rows of A
int n, //the number of columns of A
complex alpha, 
const complex x[], //the vector x
int incx, //the spacing which the elements of x occur
const complex y[], //the vector y
int incy, //the spacing which the elements of y occur
complex a[], //the complex general matrix A
int tda //the last dimension of A
);

/** f06ppc
       Perform a rank-one update for real symmetric matrix:A=alpha*x*Transpose(x)+A
Example:
	This example perform a rank-one update for A, where matrix A is
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       3.000000    7.000000    11.000000    12.000000    
       4.000000    8.000000    12.000000    16.000000    

   
    and the vector x is
       1.000000
       1.000000
       1.000000
       1.000000

void test_f06ppc()      //test for function dsyr(f06ppc)
{
    int nColRows=4;

    matrix<double> mat_input={{1,2,3,4},
                              {2,6,7,8},
                              {3,7,11,12},
                              {4,8,12,16}};

    vector<double> vec_x={1,1,1,1};

    matrix<double> mat_output;
    mat_output=mat_input;
    
    int i,j;

    dsyr(UpperTriangle,nColRows,1,vec_x,1,mat_output,nColRows);
    
    for(i=0;i<nColRows;i++)
    	for(j=0;j<i;j++)
    		mat_output[i][j]=mat_output[j][i];

    printf("Perform a rank-one update for real symmetric matrix:A=A+X*transpose(X)\n");

    printf("A=\n");                     //output mat_inputA
    for(i=0;i<nColRows;i++)
    {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_inputX
    printf("X=\n");
    for(i=0;i<nColRows;i++)
    {
      printf("       ");
      printf("%f\n",vec_x[i]);
    }

    printf("Result matrix A=\n");       //output mat_outputA
    for(i=0;i<nColRows;i++)
    {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("%f    ",mat_output[i][j]);
       printf("\n");
    }
}


   The output is as following:
Perform a rank-one update for real symmetric matrix:A=A+X*transpose(X)
A=
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       3.000000    7.000000    11.000000    12.000000    
       4.000000    8.000000    12.000000    16.000000    

X=
       1.000000
       1.000000
       1.000000
       1.000000
Result matrix A=
       2.000000    3.000000    4.000000    5.000000    
       3.000000    7.000000    8.000000    9.000000    
       4.000000    8.000000    12.000000    13.000000    
       5.000000    9.000000    13.000000    17.000000    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
Return: NULL
*/
void  dsyr(
MatrixTriangle uplo, 
int n, //the order of A
double alpha, 
const double x[], //vector x
int incx, //the spacing which the elements of x occur
double a[], //Input: the real symmetric matrix A; Output: matrix A after update
int tda //the last dimension of A
);

/** f06spc
       Perform a rank-one update for complex Hermitian matrix:A=alpha*x*ConjugateTranspose(x)+A
Example:
	This example perform a rank-one update for A, where matrix A is
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (8.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    
   
    and the vector x is
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)

void test_f06spc()    //test for function zher(f06spc)
{
   int nColRows=4;

   matrix<complex>  mat_input={{1,8+3i,2+3i,3+4i},
                               {8-3i,3,5+2i,4+4i},
                               {8-3i,5-2i,2,5+4i},
                               {3-4i,4-4i,5-4i,2}};
                               

   vector<complex>  vec_x={2+2i,1+1i,3+2i,2+4i};

   matrix<complex>  mat_output();
   mat_output=mat_input;
   
   int i,j;

   zher(UpperTriangle,nColRows,1,vec_x,1,mat_output,nColRows);
   
   for(i=0;i<nColRows;i++)
    	for(j=0;j<i;j++)
    	{
    		mat_output[i][j].m_re=mat_output[j][i].m_re;
    		mat_output[i][j].m_im=-mat_output[j][i].m_im;
    	}

   printf("Perform a rank-one update for complex Hermitian matrix:A=A+X*transpose(X)\n");

   printf("A=\n");                     //output mat_inputA
   for(i=0;i<nColRows;i++)
   {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_inputX
   printf("X=\n");
   for(i=0;i<nColRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_x[i].m_re,vec_x[i].m_im);
   }

   printf("Result matrix A=\n");                     //output mat_outputA
   for(i=0;i<nColRows;i++)
   {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("(%f,%fi)    ",mat_output[i][j].m_re,mat_output[i][j].m_im);
       printf("\n");
   }
}


   The output is as following:
Perform a rank-one update for complex Hermitian matrix:A=A+X*transpose(X)
A=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (8.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    

X=
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)
Result matrix A=
       (9.000000,0.000000i)    (12.000000,3.000000i)    (12.000000,5.000000i)    (15.000000,0.000000i)    
       (12.000000,-3.000000i)    (5.000000,0.000000i)    (10.000000,3.000000i)    (10.000000,2.000000i)    
       (12.000000,-5.000000i)    (10.000000,-3.000000i)    (15.000000,0.000000i)    (19.000000,-4.000000i)    
       (15.000000,0.000000i)    (10.000000,-2.000000i)    (19.000000,4.000000i)    (22.000000,0.000000i)    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
Return: NULL
*/
void  zher(
MatrixTriangle uplo, 
int n, //the order of A
double alpha, 
const complex x[], //the vector x 
int incx, //the spacing which the elements of x occur
complex a[], //Input: the complex Hermitian matrix A; Output: matrix A after update
int tda  //the last dimension of A
);

/** f06pqc
       Perform a rank-one update for real symmetric packed matrix:A=alpha*x*Transpose(x)+A
Example:
	This example perform a rank-one update for A, where matrix A is
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       2.000000    7.000000    9.000000    10.000000    
       4.000000    8.000000    10.000000    11.000000    

    and the vector x is
       12.000000
       13.000000
       14.000000
       15.000000

void test_f06pqc()         //test for function dspr(f06pqc)
{
    int nColRow=4;

    matrix<double> mat_input_pack={1,2,3,4,
                                     6,7,8,
                                       9,10,
                                         11};

    matrix<double> mat_input={{1,2,3,4},
                              {2,6,7,8},
                              {2,7,9,10},
                              {4,8,10,11}};

    vector<double> vec_input={12,13,14,15};

    matrix<double> mat_output_pack(mat_input_pack);
    matrix<double> mat_output(nColRow,nColRow);
    

    dspr(UpperTriangle,nColRow,1,vec_input,1,mat_output_pack);
    
    int i,j;
    int k=0;
    for(i=0;i<nColRow;i++)
    	for(j=i;j<nColRow;j++)
    		mat_output[i][j]=mat_output_pack[0][k++];
    	
    for(i=0;i<nColRow;i++)
    	for(j=0;j<i;j++)
    		mat_output[i][j]=mat_output[j][i];
    		

    printf("Perform a rank-one update for real symmetric packed matrix:A=A+X*transpose(X)\n");

    printf("A=\n");                     //output mat_input
    for(i=0;i<nColRow;i++)
    {
       printf("       ");
       for(j=0;j<nColRow;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_input
    printf("X=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_input[i]);
    }

    printf("\n");                       //output mat_output
    printf("Result matrix A=\n");
    for(i=0;i<nColRow;i++)
    {
        for(j=0;j<nColRow;j++)
        {
        	printf("       ");
            printf("%f",mat_output[i][j]);
        }
        printf("\n");
    }
}


   The output is as following:
Perform a rank-one update for real symmetric packed matrix:A=A+X*transpose(X)
A=
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       2.000000    7.000000    9.000000    10.000000    
       4.000000    8.000000    10.000000    11.000000    

X=
       12.000000
       13.000000
       14.000000
       15.000000

Result matrix A=
       145.000000       158.000000       171.000000       184.000000
       158.000000       175.000000       189.000000       203.000000
       171.000000       189.000000       205.000000       220.000000
       184.000000       203.000000       220.000000       236.000000

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
Return: NULL
*/
void  dspr(
MatrixTriangle uplo, 
int n, //the order of A
double alpha,
const double x[], //the vector x
int incx, //the spacing which the elements of x occur
double ap[] //Input: the real symmetrix packed matrix A; Output: matrix A after update
);


/** f06sqc
       Perform a rank-one update for complex Hermitian packed matrix:A=alpha*x*ConjugateTranspose(x)+A
Example:
	This example perform a rank-one update for A, where matrix A is
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (8.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    

    and the vector x is
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)


void test_f06sqc()    //test for function zhpr(f06sqc)
{
   int nColRows=4;

   matrix<complex>  mat_input_pack={1,8+3i,2+3i,3+4i,3,5+2i,4+4i,2,5+4i,2};
   
   matrix<complex>  mat_input={{1,8+3i,2+3i,3+4i},
                               {8-3i,3,5+2i,4+4i},
                               {8-3i,5-2i,2,5+4i},
                               {3-4i,4-4i,5-4i,2}};
                               

   vector<complex>  vec_x={2+2i,1+1i,3+2i,2+4i};

   matrix<complex>  mat_output_pack(mat_input_pack);
   
   
   int i,j;

   zhpr(UpperTriangle,nColRows,1,vec_x,1,mat_output_pack);
   
   int k=0;
   matrix<complex> mat_output(nColRows,nColRows);
   for(i=0;i<nColRows;i++)
   	for(j=i;j<nColRows;j++)
   	{
   		mat_output[i][j].m_re=mat_output_pack[0][k].m_re;
   		mat_output[i][j].m_im=mat_output_pack[0][k++].m_im;
   	}
   
   for(i=0;i<nColRows;i++)
   	for(j=0;j<i;j++)
   	{
   		mat_output[i][j].m_re=mat_output[j][i].m_re;
   		mat_output[i][j].m_im=-mat_output[j][i].m_im;
   	}   

   printf("Perform a rank-one update for complex Hermitian packed matrix:A=A+X*transpose(X)\n");

   printf("A=\n");                     //output mat_inputA
   for(i=0;i<nColRows;i++)
   {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_inputX
   printf("X=\n");
   for(i=0;i<nColRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_x[i].m_re,vec_x[i].m_im);
   }

   printf("Result matrix A=\n");                     //output mat_outputA
   for(i=0;i<nColRows;i++)
   {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("(%f,%fi)    ",mat_output[i][j].m_re,mat_output[i][j].m_im);
       printf("\n");
   }
}


   The output is as following:
Perform a rank-one update for complex Hermitian packed matrix:A=A+X*transpose(X)
A=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (8.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    

X=
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)
Result matrix A=
       (9.000000,0.000000i)    (12.000000,3.000000i)    (12.000000,5.000000i)    (15.000000,0.000000i)    
       (12.000000,-3.000000i)    (5.000000,0.000000i)    (10.000000,3.000000i)    (10.000000,2.000000i)    
       (12.000000,-5.000000i)    (10.000000,-3.000000i)    (15.000000,0.000000i)    (19.000000,-4.000000i)    
       (15.000000,0.000000i)    (10.000000,-2.000000i)    (19.000000,4.000000i)    (22.000000,0.000000i)    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
Return: NULL
*/
void  zhpr(
MatrixTriangle uplo, 
int n, //the order of A
double alpha,
const complex x[], //the vector x
int incx, //the spacing which the elements of x occur
complex ap[] //Input: the complex Hermitian matrix A; Output: matrix A after update
);


/** f06prc
       Perform a rank-two update for real symmetric matrix:A=alpha*x*Transpose(y)+alpha*y*Transpose(x)+A
Example:
	This example perform a rank-two update for A, where matrix A is
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       3.000000    7.000000    11.000000    12.000000    
       4.000000    8.000000    12.000000    16.000000  
       
    and the vector x is
       1.000000
       1.000000
       1.000000
       1.000000
       
    and the vector y is
       2.000000
       2.000000
       2.000000
       2.000000

void test_f06prc()      //test for function dsyr2(f06prc)
{
    int nColRows=4;

    matrix<double> mat_input={{1,2,3,4},
                              {2,6,7,8},
                              {3,7,11,12},
                              {4,8,12,16}};

    vector<double> vec_x={1,1,1,1};
    vector<double> vec_y={2,2,2,2};

    matrix<double> mat_output;
    mat_output=mat_input;
    
    int i,j;

    dsyr2(UpperTriangle,nColRows,1,vec_x,1,vec_y,1,mat_output,nColRows);
    
    
    for(i=0;i<nColRows;i++)
    	for(j=0;j<i;j++)
    		mat_output[i][j]=mat_output[j][i];


    printf("Perform a rank-two update for real symmetric matrix:A=A+X*transpose(Y)+Y*transpose(X)\n");

    printf("A=\n");                     //output mat_inputA
    for(i=0;i<nColRows;i++)
    {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_inputX
    printf("X=\n");
    for(i=0;i<nColRows;i++)
    {
      printf("       ");
      printf("%f\n",vec_x[i]);
    }
    
    printf("\n");                       //output vec_inputY
    printf("Y=\n");
    for(i=0;i<nColRows;i++)
    {
      printf("       ");
      printf("%f\n",vec_y[i]);
    }    

    printf("Result matrix A=\n");       //output mat_outputA
    for(i=0;i<nColRows;i++)
    {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("%f    ",mat_output[i][j]);
       printf("\n");
    }
}

   The output is as following:
Perform a rank-two update for real symmetric matrix:A=A+X*transpose(Y)+Y*transpose(X)
A=
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       3.000000    7.000000    11.000000    12.000000    
       4.000000    8.000000    12.000000    16.000000    

X=
       1.000000
       1.000000
       1.000000
       1.000000

Y=
       2.000000
       2.000000
       2.000000
       2.000000
Result matrix A=
       5.000000    6.000000    7.000000    8.000000    
       6.000000    10.000000    11.000000    12.000000    
       7.000000    11.000000    15.000000    16.000000    
       8.000000    12.000000    16.000000    20.000000    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
Return: NULL
*/
void  dsyr2(
MatrixTriangle uplo, 
int n, //the order of A
double alpha, 
const double x[], //the vector x
int incx, //the spacing which the elements of x occur
const double y[], //the vector y
int incy, //the spacing which the elements of y occur
double a[], //Input: the real symmetric matrix A; Output: matrix A after update
int tda //the last dimension of A
);

/** f06src
       Perform a rank-two update for complex Hermitian matrix:A=alpha*x*ConjugateTranspose(y)+alpha*y*ConjugateTranspose(x)+A
Example:
	This example perform a rank-two update for A, where matrix A is
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (8.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    

       
    and the vector x is
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)

       
    and the vector y is
       (1.000000,1.000000i)
       (2.000000,2.000000i)
       (3.000000,1.000000i)
       (4.000000,1.000000i)


void test_f06src()    //test for function zher2(f06src)
{
   int nColRows=4;

   matrix<complex>  mat_input={{1,8+3i,2+3i,3+4i},
                               {8-3i,3,5+2i,4+4i},
                               {8-3i,5-2i,2,5+4i},
                               {3-4i,4-4i,5-4i,2}};
                               

   vector<complex>  vec_x={2+2i,1+1i,3+2i,2+4i};
   vector<complex>  vec_y={1+1i,2+2i,3+1i,4+1i};

   matrix<complex>  mat_output();
   mat_output=mat_input;
   
   int i,j;

   zher2(UpperTriangle,nColRows,1,vec_x,1,vec_y,1,mat_output,nColRows);
   
    
   for(i=0;i<nColRows;i++)
    	for(j=0;j<i;j++)
    	{
    		mat_output[i][j].m_re=mat_output[j][i].m_re;
    		mat_output[i][j].m_im=-mat_output[j][i].m_im;
    	}


   printf("Perform a rank-two update for complex Hermitian matrix:A=A+X*transpose(Y)+Y*transpose(X)\n");

   printf("A=\n");                     //output mat_inputA
   for(i=0;i<nColRows;i++)
   {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_inputX
   printf("X=\n");
   for(i=0;i<nColRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_x[i].m_re,vec_x[i].m_im);
   }

   printf("\n");                       //output vec_inputY
   printf("Y=\n");
   for(i=0;i<nColRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_y[i].m_re,vec_y[i].m_im);
   }   
   
   printf("Result matrix A=\n");                     //output mat_outputA
   for(i=0;i<nColRows;i++)
   {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("(%f,%fi)    ",mat_output[i][j].m_re,mat_output[i][j].m_im);
       printf("\n");
   }
}

   The output is as following:
Perform a rank-two update for complex Hermitian matrix:A=A+X*transpose(Y)+Y*transpose(X)
A=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (8.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    

X=
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)

Y=
       (1.000000,1.000000i)
       (2.000000,2.000000i)
       (3.000000,1.000000i)
       (4.000000,1.000000i)
Result matrix A=
       (9.000000,0.000000i)    (18.000000,3.000000i)    (15.000000,8.000000i)    (19.000000,8.000000i)    
       (18.000000,-3.000000i)    (11.000000,0.000000i)    (19.000000,6.000000i)    (21.000000,3.000000i)    
       (15.000000,-8.000000i)    (19.000000,-6.000000i)    (24.000000,0.000000i)    (29.000000,-1.000000i)    
       (19.000000,-8.000000i)    (21.000000,-3.000000i)    (29.000000,1.000000i)    (26.000000,0.000000i)    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
Return: NULL
*/

void  zher2(
MatrixTriangle uplo, 
int n, //the order of A
complex alpha, 
const complex x[], //the vector x
int incx, //the spacing which the elements of x occur
const complex y[],//the vector y
int incy, //the spacing which the elements of y occur
complex a[], //Input: the complex Hermitian matrix A; Output: matrix A after update  
int tda //the last dimension of A
);

/** f06psc
       Perform a rank-two update for real symmetric packed matrix:A=alpha*x*Transpose(y)+alpha*y*Transpose(x)+A
Example:
	This example perform a rank-two update for A, where matrix A is
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       2.000000    7.000000    9.000000    10.000000    
       4.000000    8.000000    10.000000    11.000000  
       
    and the vector x is
       12.000000
       13.000000
       14.000000
       15.000000
       
    and the vector y is
       1.000000
       1.000000
       1.000000
       1.000000

void test_f06psc()         //test for function dspr2(f06psc)
{
    int nColRow=4;

    matrix<double> mat_input_pack={1,2,3,4,
                                     6,7,8,
                                       9,10,
                                         11};

    matrix<double> mat_input={{1,2,3,4},
                              {2,6,7,8},
                              {2,7,9,10},
                              {4,8,10,11}};

    vector<double> vec_x={12,13,14,15};
    vector<double> vec_y={1,1,1,1};

    matrix<double> mat_output_pack(mat_input_pack);
    matrix<double> mat_output(nColRow,nColRow);
    

    dspr2(UpperTriangle,nColRow,1,vec_x,1,vec_y,1,mat_output_pack);
    
    int i,j;
    int k=0;
    for(i=0;i<nColRow;i++)
    	for(j=i;j<nColRow;j++)
    		mat_output[i][j]=mat_output_pack[0][k++];
    	
    for(i=0;i<nColRow;i++)
    	for(j=0;j<i;j++)
    		mat_output[i][j]=mat_output[j][i];
    		

    printf("Perform a rank-two update for real symmetric packed matrix:A=A+X*transpose(Y)+Y*transpose(X)\n");

    printf("A=\n");                     //output mat_input
    for(i=0;i<nColRow;i++)
    {
       printf("       ");
       for(j=0;j<nColRow;j++)
          printf("%f    ",mat_input[i][j]);
       printf("\n");
    }

    printf("\n");                       //output vec_X
    printf("X=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_x[i]);
    }
    
    printf("\n");                       //output vec_Y
    printf("Y=\n");
    for(i=0;i<nColRow;i++)
    {
      printf("       ");
      printf("%f\n",vec_y[i]);
    }

    printf("\n");                       //output mat_output
    printf("Result matrix A=\n");
    for(i=0;i<nColRow;i++)
    {
        for(j=0;j<nColRow;j++)
        {
        	printf("       ");
            printf("%f",mat_output[i][j]);
        }
        printf("\n");
    }
}

   The output is as following:
Perform a rank-two update for real symmetric packed matrix:A=A+X*transpose(Y)+Y*transpose(X)
A=
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       2.000000    7.000000    9.000000    10.000000    
       4.000000    8.000000    10.000000    11.000000    

X=
       12.000000
       13.000000
       14.000000
       15.000000

Y=
       1.000000
       1.000000
       1.000000
       1.000000

Result matrix A=
       25.000000       27.000000       29.000000       31.000000
       27.000000       32.000000       34.000000       36.000000
       29.000000       34.000000       37.000000       39.000000
       31.000000       36.000000       39.000000       41.000000

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
Return: NULL
*/
void  dspr2(
MatrixTriangle uplo, 
int n, //the order of A
double alpha,
const double x[], //the vector x
int incx, //the spacing which the elements of x occur
const double y[], //the vector y 
int incy, //the spacing which the elements of y occur
double ap[] //Input: the real symmetric packed matrix A; Output: matrix A after update
);

/** f06ssc
       Perform a rank-two update for complex Hermitian packed matrix:A=alpha*x*ConjugateTranspose(y)+conjugate(alpha)*y*ConjugateTranspose(x)+A
Example:
	This example perform a rank-two update for A, where matrix A is
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (8.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    
       
    and the vector x is
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)
       
    and the vector y is
       (4.000000,1.000000i)
       (2.000000,1.000000i)
       (1.000000,1.000000i)
       (1.000000,1.000000i)

void test_f06ssc()    //test for function zhpr2(f06ssc)
{
   int nColRows=4;

   matrix<complex>  mat_input_pack={1,8+3i,2+3i,3+4i,3,5+2i,4+4i,2,5+4i,2};
   
   matrix<complex>  mat_input={{1,8+3i,2+3i,3+4i},
                               {8-3i,3,5+2i,4+4i},
                               {8-3i,5-2i,2,5+4i},
                               {3-4i,4-4i,5-4i,2}};
                               

   vector<complex>  vec_x={2+2i,1+1i,3+2i,2+4i};
   vector<complex>  vec_y={4+1i,2+1i,1+1i,1+1i};

   matrix<complex>  mat_output_pack(mat_input_pack);
   
   
   int i,j;

   zhpr2(UpperTriangle,nColRows,1,vec_x,1,vec_y,1,mat_output_pack);
   
   int k=0;
   matrix<complex> mat_output(nColRows,nColRows);
   for(i=0;i<nColRows;i++)
   	for(j=i;j<nColRows;j++)
   	{
   		mat_output[i][j].m_re=mat_output_pack[0][k].m_re;
   		mat_output[i][j].m_im=mat_output_pack[0][k++].m_im;
   	}
   
   for(i=0;i<nColRows;i++)
   	for(j=0;j<i;j++)
   	{
   		mat_output[i][j].m_re=mat_output[j][i].m_re;
   		mat_output[i][j].m_im=-mat_output[j][i].m_im;
   	}   

   printf("Perform a rank-one update for complex Hermitian packed matrix:A=A+X*transpose(Y)+Y*transpose(X)\n");

   printf("A=\n");                     //output mat_inputA
   for(i=0;i<nColRows;i++)
   {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("(%f,%fi)    ",mat_input[i][j].m_re,mat_input[i][j].m_im);
       printf("\n");
   }

   printf("\n");                       //output vec_inputX
   printf("X=\n");
   for(i=0;i<nColRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_x[i].m_re,vec_x[i].m_im);
   }
   
   printf("\n");                       //output vec_inputY
   printf("Y=\n");
   for(i=0;i<nColRows;i++)
   {
      printf("       ");
      printf("(%f,%fi)\n",vec_y[i].m_re,vec_y[i].m_im);
   }

   printf("Result matrix A=\n");                     //output mat_outputA
   for(i=0;i<nColRows;i++)
   {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("(%f,%fi)    ",mat_output[i][j].m_re,mat_output[i][j].m_im);
       printf("\n");
   }
}

   The output is as following:
Perform a rank-one update for complex Hermitian packed matrix:A=A+X*transpose(Y)+Y*transpose(X)
A=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (8.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    

X=
       (2.000000,2.000000i)
       (1.000000,1.000000i)
       (3.000000,2.000000i)
       (2.000000,4.000000i)

Y=
       (4.000000,1.000000i)
       (2.000000,1.000000i)
       (1.000000,1.000000i)
       (1.000000,1.000000i)
Result matrix A=
       (21.000000,0.000000i)    (19.000000,2.000000i)    (20.000000,-2.000000i)    (19.000000,-10.000000i)    
       (19.000000,-2.000000i)    (9.000000,0.000000i)    (15.000000,1.000000i)    (14.000000,-2.000000i)    
       (20.000000,2.000000i)    (15.000000,-1.000000i)    (12.000000,0.000000i)    (16.000000,1.000000i)    
       (19.000000,10.000000i)    (14.000000,2.000000i)    (16.000000,-1.000000i)    (14.000000,0.000000i)    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle is being reference.	
Return: NULL
*/
void  zhpr2(
MatrixTriangle uplo, 
int n, //the order of A
complex alpha,
const complex x[], //the vector x
int incx, //the spacing which the elements of x occur
const complex y[], //the vector y
int incy, //the spacing which the elements of y occur
complex ap[] //Input: the complex Hermitian matrix A; Output: matrix A after update
);

/** f06yac
       Compute a matrix-matrix product for two real rectangular matrices
Example:
	This example computes the product of A and B, where matrix A is
       1.000000    2.000000    3.000000    4.000000    
       5.000000    6.000000    7.000000    8.000000    
       9.000000    10.000000    11.000000    12.000000    
       13.000000    14.000000    15.000000    16.000000   
       
    and matrix B is
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    

       
void test_f06yac()      //test for function dgemm(f06yac)
{
    int nCols_b=4;
    int nColRows=4;
    int nRows_a=4;
    
    

    matrix<double> mat_a={{1,2,3,4},{5,6,7,8},{9,10,11,12},{13,14,15,16}};
    matrix<double> mat_b={{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}};
    matrix<double> mat_c(nRows_a,nCols_b);
    mat_c=0;

    dgemm(NoTranspose,NoTranspose,nRows_a,nCols_b,nColRows,1,
                mat_a,nColRows,mat_b,nCols_b,1,mat_c,nCols_b);

    printf("Compute a matrix-matrix product for general real matrix C=A*B:\n");

    int i,j;
    printf("A=\n");                     //output mat_a
    for(i=0;i<nRows_a;i++)
    {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("%f    ",mat_a[i][j]);
       printf("\n");
    }

    printf("\nB=\n");                       //output mat_b
    for(i=0;i<nColRows;i++)
    {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("%f    ",mat_b[i][j]);
       printf("\n");
    }

    printf("\nC=\n");                       //output mat_c
    for(i=0;i<nRows_a;i++)
    {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("%f    ",mat_c[i][j]);
       printf("\n");
    }
}

   The output is as following:
Compute a matrix-matrix product for general real matrix C=A*B:
A=
       1.000000    2.000000    3.000000    4.000000    
       5.000000    6.000000    7.000000    8.000000    
       9.000000    10.000000    11.000000    12.000000    
       13.000000    14.000000    15.000000    16.000000    

B=
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    

C=
       10.000000    10.000000    10.000000    10.000000    
       26.000000    26.000000    26.000000    26.000000    
       42.000000    42.000000    42.000000    42.000000    
       58.000000    58.000000    58.000000    58.000000    

Parameters:
Return: NULL
*/
void  dgemm(
MatrixTranspose transa, 
MatrixTranspose transb, 
int m, //the number of rows of C
int n, //the number of columns of C
int k, //the number of columns(rows) of A which is equal to the number of rows(columns) of B
double alpha,
const double a[], //the m by k real general matrix A(with transa=Notranspose) or the k by m matrix A (with transa=Transpose)
int tda, //the last dimension of A
const double b[],//the k by n real general matrix B(with transb=Notranspose) or the n by k matrix B(with transb=Transpose)
int tdb, //the last dimension of B
double beta, 
double c[], //Input: the m by n matrix C; Output: C=beta*C+alpha*transa(A)*transb(B)
int tdc //the last dimension of C
);

/** f06zac
       Compute a matrix-matrix product for two complex rectangular matrices
Example:
	This example computes the product of A and B, where matrix A is
       (1.000000,2.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (2.000000,1.000000i)    (3.000000,1.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (1.000000,3.000000i)    (2.000000,1.000000i)    (2.000000,3.000000i)    (5.000000,4.000000i)    
       (1.000000,4.000000i)    (2.000000,3.000000i)    (6.000000,1.000000i)    (2.000000,4.000000i)    
       
    and matrix B is
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    


       
void test_f06zac()    //test for function zgemm(f06zac)
{
   int nRows_a=4;
   int nCols_b=4;
   int nColRows=4;

   matrix<complex>  mat_a={{1+2i,8+3i,2+3i,3+4i},
                           {2+1i,3+1i,5+2i,4+4i},
                           {1+3i,2+1i,2+3i,5+4i},
                           {1+4i,2+3i,6+1i,2+4i}};
                           
   matrix<complex> mat_b={{1+2i,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}}; 
   
   matrix<complex> mat_c(nRows_a,nColRows);
   mat_c=0;

   zgemm(NoTranspose,NoTranspose,nRows_a,nCols_b,nColRows,1,
                mat_a,nColRows,mat_b,nCols_b,1,mat_c,nCols_b);

   printf("Compute a matrix-matrix product for general complex matrix C=A*B:\n");

   int i,j;
   
   printf("A=\n");                     //output mat_a
   for(i=0;i<nRows_a;i++)
   {
       printf("       ");
       for(j=0;j<nColRows;j++)
          printf("(%f,%fi)    ",mat_a[i][j].m_re,mat_a[i][j].m_im);
       printf("\n");
   }

   printf("B=\n");                     //output mat_b
   for(i=0;i<nColRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_b[i][j].m_re,mat_b[i][j].m_im);
       printf("\n");
   }
   
   printf("C=\n");                     //output mat_c
   for(i=0;i<nRows_a;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_c[i][j].m_re,mat_c[i][j].m_im);
       printf("\n");
   }   
}

   The output is as following:
Compute a matrix-matrix product for general complex matrix C=A*B:
A=
       (1.000000,2.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (2.000000,1.000000i)    (3.000000,1.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (1.000000,3.000000i)    (2.000000,1.000000i)    (2.000000,3.000000i)    (5.000000,4.000000i)    
       (1.000000,4.000000i)    (2.000000,3.000000i)    (6.000000,1.000000i)    (2.000000,4.000000i)    
B=
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
C=
       (10.000000,14.000000i)    (14.000000,12.000000i)    (14.000000,12.000000i)    (14.000000,12.000000i)    
       (12.000000,12.000000i)    (14.000000,8.000000i)    (14.000000,8.000000i)    (14.000000,8.000000i)    
       (4.000000,13.000000i)    (10.000000,11.000000i)    (10.000000,11.000000i)    (10.000000,11.000000i)    
       (3.000000,14.000000i)    (11.000000,12.000000i)    (11.000000,12.000000i)    (11.000000,12.000000i)    


Parameters:
Return: NULL
*/
void  zgemm(
MatrixTranspose transa, 
MatrixTranspose transb, 
int m, //the number of rows of C
int n, //the number of columns of C
int k, //the number of columns(rows) of A, which is equal to the number of rows(columns) of B
complex alpha,
const complex a[], //the m by k complex matrix A(with transa=Notranspose) or k by m complex matrix A(with transa=Transpose or ConjugateTranspose)
int tda, //the last dimension of A
const complex b[],//the k by n complex matrix B(with transb=Notranspose) or n by k complex matrix B (with transb=Transpose or ConjugateTranspose)
int tdb, //the last dimension of B
complex beta, 
complex c[],//Input: the m by n complex matrix C; Output: C=beta*C+alpha*transa(A)*transb(B)
int tdc //the last dimension of C
);


/** f06ycc
       Compute a matrix-matrix product for one real symmetric rectangular matrix A and one real general matrix B
Example:
	This example computes the product of A and B, where the symmetric matrix A is
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       3.000000    7.000000    11.000000    12.000000    
       4.000000    8.000000    12.000000    16.000000  
       
    and the general matrix B is
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       
void test_f06ycc()      //test for function dsymm(f06ycc)
{
    int nCols_b=4;
    int nColRows_sy=4;
    
    matrix<double> mat_sy={{1,2,3,4},{2,6,7,8},{3,7,11,12},{4,8,12,16}};
    matrix<double> mat_g={{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}};
    matrix<double> mat_c(nColRows_sy,nCols_b);
    mat_c=0;

    dsymm(LeftSide,UpperTriangle,nColRows_sy,nCols_b,1,
                mat_sy,nColRows_sy,mat_g,nCols_b,1,mat_c,nCols_b);

    printf("Compute a matrix-matrix product for a symmetric real matrix A and a general real matrix B, C=A*B:\n");

    int i,j;
    printf("A=\n");                     //output symmetric real matrix mat_sy
    for(i=0;i<nColRows_sy;i++)
    {
       printf("       ");
       for(j=0;j<nColRows_sy;j++)
          printf("%f    ",mat_sy[i][j]);
       printf("\n");
    }

    printf("\nB=\n");                       //output general real mat_g
    for(i=0;i<nColRows_sy;i++)
    {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("%f    ",mat_g[i][j]);
       printf("\n");
    }

    printf("\nC=\n");                       //output mat_c
    for(i=0;i<nColRows_sy;i++)
    {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("%f    ",mat_c[i][j]);
       printf("\n");
    }
}

   The output is as following:
Compute a matrix-matrix product for a symmetric real matrix A and a general real matrix B, C=A*B:
A=
       1.000000    2.000000    3.000000    4.000000    
       2.000000    6.000000    7.000000    8.000000    
       3.000000    7.000000    11.000000    12.000000    
       4.000000    8.000000    12.000000    16.000000    

B=
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    

C=
       10.000000    10.000000    10.000000    10.000000    
       23.000000    23.000000    23.000000    23.000000    
       33.000000    33.000000    33.000000    33.000000    
       40.000000    40.000000    40.000000    40.000000    

Parameters:
	side=Leftside: compute A*B 
	side=Rightside: compute B*A 
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of A is being reference.	

Return: NULL
*/
void  dsymm(
OperationSide side, 
MatrixTriangle uplo, 
int m, //the number of rows of C
int n, //the number of columns of C
double alpha,
const double a[], //the order m symmetric matrix A(with side=Leftside) or the order n symmetric matrix A (with side=Rightside)
int tda, //the last dimension of A
const double b[], //the m by n matrix B(with side=Leftside) or the n by m matrix B(with side=Rightside)
int tdb, //the last dimension of B
double beta, 
double c[], //Input: the m by n matrix C; Output: the product of A and B
int tdc  //the last dimension of C
);

/** f06zcc
       Compute a matrix-matrix product for one complex Hermitian matrix A and one complex general matrix B
Example:
	This example computes the product of A and B, where the complex Hermitian matrix A is
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (2.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    

       
    and the general complex matrix B is
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    

       
void test_f06zcc()    //test for function zhemm(f06zcc)
{
   int nColRows_h=4;
   int nCols_b=4;

   matrix<complex>  mat_h={{1,8+3i,2+3i,3+4i},
                           {8-3i,3,5+2i,4+4i},
                           {2-3i,5-2i,2,5+4i},
                           {3-4i,4-4i,5-4i,2}};
                           
   matrix<complex> mat_b={{1+2i,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}}; 
   
   matrix<complex> mat_c(nColRows_h,nCols_b);
   mat_c=0;

   zhemm(LeftSide,UpperTriangle,nColRows_h,nCols_b,1,
                mat_h,nColRows_h,mat_b,nCols_b,1,mat_c,nCols_b);

   printf("Compute a matrix-matrix product for a complex Hermitian matrix H and a general complex matrix B, C=H*B:\n");

   int i,j;
   
   printf("H=\n");                     //output complex Hermitian matrix mat_h
   for(i=0;i<nColRows_h;i++)
   {
       printf("       ");
       for(j=0;j<nColRows_h;j++)
          printf("(%f,%fi)    ",mat_h[i][j].m_re,mat_h[i][j].m_im);
       printf("\n");
   }

   printf("B=\n");                     //output general complex matrix mat_b
   for(i=0;i<nColRows_h;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_b[i][j].m_re,mat_b[i][j].m_im);
       printf("\n");
   }
   
   printf("C=\n");                     //output mat_c
   for(i=0;i<nColRows_h;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_c[i][j].m_re,mat_c[i][j].m_im);
       printf("\n");
   }   
}

   The output is as following:
Compute a matrix-matrix product for a complex Hermitian matrix H and a general complex matrix B, C=H*B:
H=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (2.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    
B=
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
C=
       (14.000000,12.000000i)    (14.000000,10.000000i)    (14.000000,10.000000i)    (14.000000,10.000000i)    
       (26.000000,19.000000i)    (20.000000,3.000000i)    (20.000000,3.000000i)    (20.000000,3.000000i)    
       (20.000000,3.000000i)    (14.000000,-1.000000i)    (14.000000,-1.000000i)    (14.000000,-1.000000i)    
       (22.000000,-6.000000i)    (14.000000,-12.000000i)    (14.000000,-12.000000i)    (14.000000,-12.000000i)    

Parameters:
	side=Leftside: compute A*B 
	side=Rightside: compute B*A 
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of A is being reference.	

Return: NULL
*/
void  zhemm(
OperationSide side, 
MatrixTriangle uplo, 
int m, //the number of rows of C
int n, //the number of column of C
complex alpha,
const complex a[], //the complex Hermitian matrix A of m by m (with side=Leftside) or n by n (with side=Rightside)
int tda, //the last dimension of A
const complex b[], //the complex matrix B of m by n (with side=Leftside) or n by m (with side=Rightside)
int tdb, //the last dimension of B
complex beta, 
complex c[], //Input: the m by n matrix; Output: C=beta*C+alpha*AB
int tdc //the last dimension of C
);

/**	f06yfc
       Compute a matrix-matrix product for one real triangular matrix A and one real general matrix B
Example:
	This example computes the product of A and B, where the real triangular matrix A is
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    11.000000    12.000000    
       0.000000    0.000000    0.000000    16.000000    

       
    and the general real matrix B is
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       
void test_f06yfc()      //test for function dtrmm(f06yfc)
{
    int nCols_g=4;
    int nColRows_tr=4;
    
    matrix<double> mat_tr={{1,2,3,4},{0,6,7,8},{0,0,11,12},{0,0,0,16}};
    matrix<double> mat_g={{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}};
    matrix<double> mat_c(nColRows_tr,nCols_g);
    mat_c=mat_g;

    dtrmm(LeftSide,UpperTriangle,NoTranspose,NotUnitTriangular,nColRows_tr,nCols_g,1,
                mat_tr,nColRows_tr,mat_c,nCols_g);

    printf("Compute a matrix-matrix product for a triangular real matrix A and a general real matrix B, C=A*B:\n");

    int i,j;
    printf("A=\n");                     //output triangular real matrix mat_tr
    for(i=0;i<nColRows_tr;i++)
    {
       printf("       ");
       for(j=0;j<nColRows_tr;j++)
          printf("%f    ",mat_tr[i][j]);
       printf("\n");
    }

    printf("\n");
    printf("B=\n");                //output general real mat_g
    for(i=0;i<nColRows_tr;i++)
    {
       printf("       ");
       for(j=0;j<nCols_g;j++)
          printf("%f    ",mat_g[i][j]);
       printf("\n");
    }

    printf("\n");
    printf("C=\n");                    //output mat_c
    for(i=0;i<nColRows_tr;i++)
    {
       printf("       ");
       for(j=0;j<nCols_g;j++)
          printf("%f    ",mat_c[i][j]);
       printf("\n");
    }
}

   The output is as following:
Compute a matrix-matrix product for a triangular real matrix A and a general real matrix B, C=A*B:
A=
       1.000000    2.000000    3.000000    4.000000    
       0.000000    6.000000    7.000000    8.000000    
       0.000000    0.000000    11.000000    12.000000    
       0.000000    0.000000    0.000000    16.000000    

B=
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    

C=
       10.000000    10.000000    10.000000    10.000000    
       21.000000    21.000000    21.000000    21.000000    
       23.000000    23.000000    23.000000    23.000000    
       16.000000    16.000000    16.000000    16.000000    

Parameters:
	side=Leftside: compute transa(A)*B 
	side=Rightside: compute B*transa(A) 
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of A is being reference.	
    diag specify whether or not the matrix A is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.
Return: NULL
*/
void  dtrmm(
MatrixTriangle side, 
MatrixTriangle uplo, 
MatrixTranspose transa, 
MatrixUnitTriangular diag,
int m, //the number of rows of B
int n, //the number of columns of B
double alpha,
const double a[], //the real triangular matrix A of m by m(with side=Leftside) or n by n(with side=Rightside)
int tda, //the last dimension of A
double b[], //Input: the m by n general real matrix B; Output: alpha*AB
int tdb  //the last dimension of B
);

/** f06zfc
       Compute a matrix-matrix product for one complex triangular matrix A and one complex general matrix B
Example:
	This example computes the product of T and B, where the complex triangular matrix T is
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (0.000000,0.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (2.000000,0.000000i)    
  

       
    and the general complex matrix B is
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
 
       
void test_f06zfc()    //test for function ztrmm(f06zfc)
{
   int nColRows_tr=4;
   int nCols_b=4;

   matrix<complex>  mat_tr={{1,8+3i,2+3i,3+4i},
                           {0,3,5+2i,4+4i},
                           {0,0,2,5+4i},
                           {0,0,0,2}};
                           
   matrix<complex> mat_b={{1+2i,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}}; 
   
   matrix<complex> mat_c(nColRows_tr,nCols_b);
   mat_c=mat_b;

   ztrmm(LeftSide,UpperTriangle,NoTranspose,NotUnitTriangular,nColRows_tr,nCols_b,1,
                mat_tr,nColRows_tr,mat_c,nCols_b);

   printf("Compute a matrix-matrix product for a complex triangular matrix H and a general complex matrix B, C=T*B:\n");

   int i,j;
   
   printf("T=\n");                     //output complex triangular matrix mat_tr
   for(i=0;i<nColRows_tr;i++)
   {
       printf("       ");
       for(j=0;j<nColRows_tr;j++)
          printf("(%f,%fi)    ",mat_tr[i][j].m_re,mat_tr[i][j].m_im);
       printf("\n");
   }

   printf("B=\n");                     //output general complex matrix mat_b
   for(i=0;i<nColRows_tr;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_b[i][j].m_re,mat_b[i][j].m_im);
       printf("\n");
   }
   
   printf("C=\n");                     //output mat_c
   for(i=0;i<nColRows_tr;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_c[i][j].m_re,mat_c[i][j].m_im);
       printf("\n");
   }   
}

   The output is as following:
Compute a matrix-matrix product for a complex triangular matrix H and a general complex matrix B, C=T*B:
T=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (0.000000,0.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (2.000000,0.000000i)    
B=
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
C=
       (14.000000,12.000000i)    (14.000000,10.000000i)    (14.000000,10.000000i)    (14.000000,10.000000i)    
       (12.000000,6.000000i)    (12.000000,6.000000i)    (12.000000,6.000000i)    (12.000000,6.000000i)    
       (7.000000,4.000000i)    (7.000000,4.000000i)    (7.000000,4.000000i)    (7.000000,4.000000i)    
       (2.000000,0.000000i)    (2.000000,0.000000i)    (2.000000,0.000000i)    (2.000000,0.000000i)    


Parameters:
	side=Leftside: compute transa(A)*B 
	side=Rightside: compute B*transa(A) 
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of A is being reference.	
    diag specify whether or not the matrix A is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/
void  ztrmm(
MatrixTriangle side, 
MatrixTriangle uplo, 
MatrixTranspose transa, 
MatrixUnitTriangular diag,
int m, //the number of rows of B
int n, //the number of columns of B
complex alpha,
const complex a[], //the complex triangular matrix A of m by m(with side=Leftside) or n by n(with side= Rightside)
int tda, //the last dimension of A
complex b[],//Input: the m by n complex general matrix B; Output: alpha*AB
int tdb //the last dimension of B
);

/** f06zjc
       Solve a system of equations with multiple right-hand sides for real triangular coefficient matrix
Example:
	This example solve a system of equations with multiple fight-hand sides, where the coefficient matrix T is
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (0.000000,0.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (2.000000,0.000000i)    
       
    and the multiple right-hand sides matrix B is
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (2.000000,3.000000i)    (3.000000,2.000000i)    (4.000000,0.000000i)    
       (2.000000,-2.000000i)    (3.000000,5.000000i)    (4.000000,0.000000i)    (5.000000,0.000000i)    
       (3.000000,-1.000000i)    (4.000000,2.000000i)    (5.000000,-3.000000i)    (6.000000,4.000000i)    
      
void test_f06zjc()    //test for function ztrsm(f06zjc)
{
   int nColRows_tr=4;
   int nCols_b=4;

   matrix<complex>  mat_tr={{1,8+3i,2+3i,3+4i},
                           {0,3,5+2i,4+4i},
                           {0,0,2,5+4i},
                           {0,0,0,2}};
                           
   matrix<complex> mat_b={{1+2i,1,1,1},{1,2+3i,3+2i,4},{2-2i,3+5i,4,5},{3-1i,4+2i,5-3i,6+4i}}; 
   
   matrix<complex> mat_c(nColRows_tr,nCols_b);
   mat_c=mat_b;

   ztrsm(LeftSide,UpperTriangle,NoTranspose,NotUnitTriangular,nColRows_tr,nCols_b,1,
                mat_tr,nColRows_tr,mat_c,nCols_b);

   printf("Solve a systrm of equations with multiple right-hand sides, T*X=B:\n");

   int i,j;
   
   printf("T=\n");                     //output complex triangular matrix mat_tr
   for(i=0;i<nColRows_tr;i++)
   {
       printf("       ");
       for(j=0;j<nColRows_tr;j++)
          printf("(%f,%fi)    ",mat_tr[i][j].m_re,mat_tr[i][j].m_im);
       printf("\n");
   }

   printf("B=\n");                     //output right-hand complex matrix mat_b
   for(i=0;i<nColRows_tr;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_b[i][j].m_re,mat_b[i][j].m_im);
       printf("\n");
   }
   
   printf("X=\n");                     //output mat_c
   for(i=0;i<nColRows_tr;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_c[i][j].m_re,mat_c[i][j].m_im);
       printf("\n");
   }   
}

   The output is as following:
Solve a systrm of equations with multiple right-hand sides, T*X=B:
T=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (0.000000,0.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (0.000000,0.000000i)    (0.000000,0.000000i)    (0.000000,0.000000i)    (2.000000,0.000000i)    
B=
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (2.000000,3.000000i)    (3.000000,2.000000i)    (4.000000,0.000000i)    
       (2.000000,-2.000000i)    (3.000000,5.000000i)    (4.000000,0.000000i)    (5.000000,0.000000i)    
       (3.000000,-1.000000i)    (4.000000,2.000000i)    (5.000000,-3.000000i)    (6.000000,4.000000i)    
X=
       (-5.666667,-38.000000i)    (10.666667,-33.333333i)    (-38.333333,-52.000000i)    (51.333333,-74.666667i)    
       (2.083333,5.750000i)    (-0.833333,4.666667i)    (6.916667,6.250000i)    (-5.666667,12.333333i)    
       (-3.750000,-2.750000i)    (-1.500000,-4.000000i)    (-7.250000,-1.250000i)    (-1.000000,-11.000000i)    
       (1.500000,-0.500000i)    (2.000000,1.000000i)    (2.500000,-1.500000i)    (3.000000,2.000000i)    

Parameters:
	side=Leftside: solve the equations trans(A)*X=alpha*B
	side=Rightside: solve the equations X*trans(A)=alpha*B
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of A is being reference.	
    diag specify whether or not the matrix A is unit triangular. if diag=UnitTriangular, the diagonal elements are not referenced.

Return: NULL
*/
void  ztrsm(
MatrixTriangle side, 
MatrixTriangle uplo, 
MatrixTranspose transa, 
MatrixUnitTriangular diag,
int m, //the number of rows of B
int n, //the number of columns of B
complex alpha,
const complex a[], //the coefficient matrix A of m by m(with side=Leftside) or n by n (with side=Rightside)
int tda, //the last dimension of A
complex b[], //Input: the m by n right-hand side matrix B; Output:the solvtion of the equations 
int tdb //the last dimension of B
);

/** f06ypc
       Perform a rank-k update of a real symmetric matrix
Example:
	This example perform a rank-k update of a real symmetric matrix C.

void test_f06ypc()      //test for function dsyrk(f06ypc)
{
    int nCols=4;
    int nRows=4;
    
    matrix<double> mat_a={{1,2,3,4},{2,3,4,5},{3,7,11,12},{4,8,12,16}};
    matrix<double> mat_c(nRows,nRows);
    mat_c=0;

    dsyrk(UpperTriangle,NoTranspose,nRows,nCols,1,
                mat_a,nCols,1,mat_c,nRows);
                
    int i,j;
    for(i=0;i<nRows;i++)
    	for(j=0;j<i;j++)
    		mat_c[i][j]=mat_c[j][i];
    
    printf("Perform a rank-k update of a real symmtric matrix C=A*transpose(A)+C:\n");

    printf("A=\n");                     //output general real matrix mat_sy
    for(i=0;i<nCols;i++)
    {
       printf("       ");
       for(j=0;j<nRows;j++)
          printf("%f    ",mat_a[i][j]);
       printf("\n");
    }

    printf("\n"); 
    printf("C=\n");                       //output real symmtric matrix mat_c
    for(i=0;i<nRows;i++)
    {
       printf("       ");
       for(j=0;j<nRows;j++)
          printf("%f    ",mat_c[i][j]);
       printf("\n");
    }
}

   The output is as following:
Perform a rank-k update of a real symmtric matrix C=A*transpose(A)+C:
A=
       1.000000    2.000000    3.000000    4.000000    
       2.000000    3.000000    4.000000    5.000000    
       3.000000    7.000000    11.000000    12.000000    
       4.000000    8.000000    12.000000    16.000000    

C=
       30.000000    40.000000    98.000000    120.000000    
       40.000000    54.000000    131.000000    160.000000    
       98.000000    131.000000    323.000000    392.000000    
       120.000000    160.000000    392.000000    480.000000    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of C is being reference.	
    trans=Notranspose: return alpha*A*transpose(A)+beta*C in array C;
    trans=Transpose or ConjugateTranspose: return alpha*transpose(A)*A+beta*C in array C;
Return: NULL
*/
void  dsyrk(
MatrixTriangle uplo, 
MatrixTranspose trans, 
int n, //the order of C
int k, //the number of rows(with tran=Transpose) or columns(tran=Notranspose) of A
double alpha,
const double a[], //the matrix A
int tda, //the last dimension of A
double beta,
double c[], //Input: the order n real symmetric matrix C; Output: matrix C after update
int tdc //the last dimension of C
);

/** f06zpc
       Perform a rank-k update of a complex Hermitian matrix
Example:
	This example perform a rank-k update of a complex Hermitian matrix C

void test_f06zpc()    //test for function zherk(f06zpc)
{
   int nRows=4;
   int nCols=4;

   matrix<complex>  mat_a={{1,8+3i,2+3i,3+4i},
                           {8-3i,3,5+2i,4+4i},
                           {2-3i,5-2i,2,5+4i},
                           {3-4i,4-4i,5-4i,2}};
                           
   matrix<complex> mat_c(nRows,nRows);
   mat_c=0;

   zherk(UpperTriangle,NoTranspose,nRows,nCols,1,
                mat_a,nCols,1,mat_c,nRows);
                
   int i,j;
   for(i=0;i<nRows;i++)
    	for(j=0;j<i;j++)
    	{
    		mat_c[i][j].m_re=mat_c[j][i].m_re;  
    		mat_c[i][j].m_im=-mat_c[j][i].m_im; 
    	}
   
   printf("Perform a rank-k update of a complex Hermitian matrix C=H*conj(H)+C:\n");

   
   
   printf("H=\n");                     //output complex  matrix mat_h
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_a[i][j].m_re,mat_a[i][j].m_im);
       printf("\n");
   }

   
   printf("C=\n");                     //output complex Hermitian matrix  mat_c
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nRows;j++)
          printf("(%f,%fi)    ",mat_c[i][j].m_re,mat_c[i][j].m_im);
       printf("\n");
   }   
}

   The output is as following:
Perform a rank-k update of a complex Hermitian matrix C=H*conj(H)+C:
H=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (2.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    
C=
       (112.000000,0.000000i)    (76.000000,27.000000i)    (71.000000,48.000000i)    (27.000000,79.000000i)    
       (76.000000,-27.000000i)    (143.000000,0.000000i)    (86.000000,32.000000i)    (73.000000,73.000000i)    
       (71.000000,-48.000000i)    (86.000000,-32.000000i)    (87.000000,0.000000i)    (66.000000,27.000000i)    
       (27.000000,-79.000000i)    (73.000000,-73.000000i)    (66.000000,-27.000000i)    (102.000000,0.000000i)    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of C is being reference.	
    trans=Notranspose: return alpha*A*transpose(A)+beta*C in array C;
    trans=ConjugateTranspose: return alpha*ConjugateTranspose(A)*A+beta*C in array C;

Return: NULL
		
*/
void  zherk(
MatrixTriangle uplo, 
MatrixTranspose trans, 
int n, //the order of C
int k, //the number of rows(with trans=Transpose/ConjugateTranspose) or columns(trans=Notranspose) of A
double alpha,
const complex a[], //the matrix A
int tda, //the last dimension of A
double beta,
complex c[], //Input: the Hermitian matrix C; Output: the matrix C after update
int tdc  //the last dimension of A
);

/** f06yrc
       Perform a rank-2k update of a real symmetric matrix
Example:
	This example perform a rank-2k update of a real symmtric matrix C

void test_f06yrc()      //test for function dsyr2k(f06yrc)
{
    int nCols=4;
    int nRows=4;
    int nColRows=4;
    
    matrix<double> mat_a={{1,2,3,4},{2,2,8,4},{3,7,11,12},{4,8,12,16}};
    matrix<double> mat_b={{1,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}};
    matrix<double> mat_c(nRows,nRows);
    mat_c=0;

    dsyr2k(UpperTriangle,NoTranspose,nRows,nCols,1,
                mat_a,nCols,mat_b,nCols,1,mat_c,nRows);
                
    int i,j;
    for(i=0;i<nRows;i++)
    	for(j=0;j<i;j++)
    		mat_c[i][j]=mat_c[j][i];

    printf("Perform a rank-2k update of a real symmetric matrix C=C+A*transpose(B)+B*transpose(A):\n");

    
    printf("A=\n");                     //output general real matrix A
    for(i=0;i<nRows;i++)
    {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("%f    ",mat_a[i][j]);
       printf("\n");
    }

    printf("B=\n");                       //output general real matrix B
    for(i=0;i<nRows;i++)
    {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("%f    ",mat_b[i][j]);
       printf("\n");
    }

    printf("C=\n");                       //output symmetric matrix C
    for(i=0;i<nRows;i++)
    {
       printf("       ");
       for(j=0;j<nRows;j++)
          printf("%f    ",mat_c[i][j]);
       printf("\n");
    }
}

   The output is as following:
Perform a rank-2k update of a real symmetric matrix C=C+A*transpose(B)+B*transpose(A):
A=
       1.000000    2.000000    3.000000    4.000000    
       2.000000    2.000000    8.000000    4.000000    
       3.000000    7.000000    11.000000    12.000000    
       4.000000    8.000000    12.000000    16.000000    
B=
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
       1.000000    1.000000    1.000000    1.000000    
C=
       20.000000    26.000000    43.000000    50.000000    
       26.000000    32.000000    49.000000    56.000000    
       43.000000    49.000000    66.000000    73.000000    
       50.000000    56.000000    73.000000    80.000000    

Parameters:
    uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of C is being reference.	
    trans=Notranspose: both A and B are n by k, return alpha*A*transpose(B)+alpha*B*transpose(A)+beta*C in array C;
    trans=Transpose or ConjugateTranspose: both A and B are k by n, return alpha*transpose(A)*B+alpha*transpose(B)*A+beta*C in array C;

Return: NULL
		
*/
void  dsyr2k(
MatrixTriangle uplo, 
MatrixTranspose trans, 
int n, //the order of C 
int k, 
double alpha,
const double a[], //the general matrix A
int tda, //the last dimesion of A
const double b[],//the general matrix B
int tdb, //the last dimension of B
double beta, 
double c[], //Input: the order n real symmetric matrix C; Output: matrix C after update
int tdc //the last dimension of C
);

/** f06zrc
       Perform a rank-2k update of a complex Hermitian matrix
Example:
	This example perform a rank-2k update of complex Hermitian matrix C

void test_f06zrc()    //test for function zher2k(f06zrc)
{
   int nCols=4;
   int nRows=4;

   matrix<complex>  mat_a={{1,8+3i,2+3i,3+4i},
                           {8-3i,3,5+2i,4+4i},
                           {2-3i,5-2i,2,5+4i},
                           {3-4i,4-4i,5-4i,2}};
                           
   matrix<complex> mat_b={{1+2i,1,1,1},{1,1-4i,1,1},{1,1+1i,1,1},{1-4i,1,1,1}}; 
   
   matrix<complex> mat_c(nRows,nRows);
   mat_c=0;

   zher2k(UpperTriangle,NoTranspose,nRows,nCols,1,
                mat_a,nCols,mat_b,nCols,1,mat_c,nRows);


   printf("Perform a rank-2k update of a complex Hermitian matrix C=C+A*transpose(B)+B*transpose(A):\n");

   int i,j;
   for(i=0;i<nRows;i++)
    	for(j=0;j<i;j++)
    	{
    		mat_c[i][j].m_re=mat_c[j][i].m_re;
    		mat_c[i][j].m_im=-mat_c[j][i].m_im;
    	}
   
   printf("A=\n");                     //output general complex matrix mat_a
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_a[i][j].m_re,mat_a[i][j].m_im);
       printf("\n");
   }

   printf("B=\n");                     //output general complex matrix mat_b
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_b[i][j].m_re,mat_b[i][j].m_im);
       printf("\n");
   }
   
   printf("C=\n");                     //output complex Hermitian matrix mat_c
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nRows;j++)
          printf("(%f,%fi)    ",mat_c[i][j].m_re,mat_c[i][j].m_im);
       printf("\n");
   }   
}

   The output is as following:
Perform a rank-2k update of a complex Hermitian matrix C=C+A*transpose(B)+B*transpose(A):
A=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (2.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    
B=
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,-4.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,1.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,-4.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
C=
       (28.000000,0.000000i)    (16.000000,55.000000i)    (25.000000,7.000000i)    (20.000000,32.000000i)    
       (16.000000,-55.000000i)    (40.000000,0.000000i)    (42.000000,-19.000000i)    (62.000000,31.000000i)    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of C is being reference.	
    trans=Notranspose: both A and B are n by k, return alpha*A*ConjugateTranspose(B)+ConjugateTranspose(alpha)*B*ConjugateTranspose(A)+beta*C in array C;
    trans=ConjugateTranspose: both A and B are k by n, return alpha*ConjugateTranspose(A)*B+ConjugateTranspose(alpha)*ConjugateTranspose(B)*A+beta*C in array C;
    
Return: NULL
		
*/
void  zher2k(
MatrixTriangle uplo, 
MatrixTranspose trans, 
int n, //the order of C
int k, 
complex alpha,
const complex a[], //the complex general matrix A
int tda, //the last dimension of A
const complex b[], //the complex general matrix B
int tdb, //the last dimension of B
double beta, 
complex c[], //Input: the order n complex Hermitian matrix C; Output: matrix C after update
int tdc //the last dimension of C
);

/** f06ztc
       Compute a matrix-matrix product for one complex symmetric matrix and one somplex rectangular matrix
Example:
	This example computes a matrix-matrix product for A and B, where the complex symmetric matrix A is
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (2.000000,3.000000i)    (5.000000,2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,4.000000i)    (4.000000,4.000000i)    (5.000000,4.000000i)    (2.000000,0.000000i)    
   
    and the matrix B is
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
    

void test_f06ztc()    //test for function zsymm(f06ztc)
{
   int nColRows_sy=4;
   int nCols_b=4;

   matrix<complex>  mat_sy={{1,8+3i,2+3i,3+4i},
                           {8+3i,3,5+2i,4+4i},
                           {2+3i,5+2i,2,5+4i},
                           {3+4i,4+4i,5+4i,2}};
                           
   matrix<complex> mat_b={{1+2i,1,1,1},{1,1,1,1},{1,1,1,1},{1,1,1,1}}; 
   
   matrix<complex> mat_c(nColRows_sy,nCols_b);
   mat_c=0;

   zsymm(LeftSide,UpperTriangle,nColRows_sy,nCols_b,1,
                mat_sy,nColRows_sy,mat_b,nCols_b,1,mat_c,nCols_b);

   printf("Compute a matrix-matrix product for a complex symmetric matrix A and a general complex matrix B, C=A*B:\n");

   int i,j;
   
   printf("A=\n");                     //output complex symmetric matrix mat_h
   for(i=0;i<nColRows_sy;i++)
   {
       printf("       ");
       for(j=0;j<nColRows_sy;j++)
          printf("(%f,%fi)    ",mat_sy[i][j].m_re,mat_sy[i][j].m_im);
       printf("\n");
   }

   printf("B=\n");                     //output general complex matrix mat_b
   for(i=0;i<nColRows_sy;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_b[i][j].m_re,mat_b[i][j].m_im);
       printf("\n");
   }
   
   printf("C=\n");                     //output mat_c
   for(i=0;i<nColRows_sy;i++)
   {
       printf("       ");
       for(j=0;j<nCols_b;j++)
          printf("(%f,%fi)    ",mat_c[i][j].m_re,mat_c[i][j].m_im);
       printf("\n");
   }   
}

   The output is as following:
Compute a matrix-matrix product for a complex symmetric matrix A and a general complex matrix B, C=A*B:
A=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (2.000000,3.000000i)    (5.000000,2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,4.000000i)    (4.000000,4.000000i)    (5.000000,4.000000i)    (2.000000,0.000000i)    
B=
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
C=
       (14.000000,12.000000i)    (14.000000,10.000000i)    (14.000000,10.000000i)    (14.000000,10.000000i)    
       (14.000000,25.000000i)    (20.000000,9.000000i)    (20.000000,9.000000i)    (20.000000,9.000000i)    
       (8.000000,13.000000i)    (14.000000,9.000000i)    (14.000000,9.000000i)    (14.000000,9.000000i)    
       (6.000000,18.000000i)    (14.000000,12.000000i)    (14.000000,12.000000i)    (14.000000,12.000000i)    

Parameters:
	side=Leftside: compute A*B 
	side=Rightside: compute B*A 
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of A is being reference.	

Return: NULL
*/
void  zsymm(
OperationSide side, 
MatrixTriangle uplo, 
int m, //the number of rows of C
int n, //the number of columns of C
complex alpha,
const complex a[], //the real symmetric matrix A of n by n(side=Rightside) or m by m(side=Leftside)
int tda, //the last dimension of A
const complex b[],//the real matrix B of m by n(side=Leftside) or n by m (side=Rightside)
int tdb, //the last dimension of B
complex beta, 
complex c[],//Input: the m by n matrix; Output: alpha*AB+beta*C
int tdc  //the last dimension of C
);

/** f06zuc
       Perform a rank-k update of a complex symmetric matrix 
Example:
	This example perform a rank-k update of a complex sysmetric matrix C
   

void test_f06zuc()    //test for function zsyrk(f06zuc)
{
   int nRows=4;
   int nCols=4;

   matrix<complex>  mat_a={{1,8+3i,2+3i,3+4i},
                           {8+3i,3,5+2i,4+4i},
                           {2+3i,5-2i,2,5+4i},
                           {3+4i,4-4i,5-4i,2}};
                           
   matrix<complex> mat_c(nRows,nRows);
   mat_c=0;

   zsyrk(UpperTriangle,NoTranspose,nRows,nCols,1,
                mat_a,nCols,1,mat_c,nRows);

   printf("Perform a rank-k update of a complex symmetric matrix C=A*conj(A)+C:\n");

   int i,j;
   for(i=0;i<nRows;i++)
    	for(j=0;j<i;j++)
    	{
    		mat_c[i][j].m_re=mat_c[j][i].m_re;
    		mat_c[i][j].m_im=mat_c[j][i].m_im;
    	}
   
   printf("A=\n");                     //output complex  matrix mat_h
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_a[i][j].m_re,mat_a[i][j].m_im);
       printf("\n");
   }

   
   printf("C=\n");                     //output complex Hermitian matrix  mat_c
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nRows;j++)
          printf("(%f,%fi)    ",mat_c[i][j].m_re,mat_c[i][j].m_im);
       printf("\n");
   }   
}

   The output is as following:
Perform a rank-k update of a complex symmetric matrix C=A*conj(A)+C:
A=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (2.000000,3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    
C=
       (44.000000,84.000000i)    (32.000000,59.000000i)    (51.000000,40.000000i)    (75.000000,-1.000000i)    
       (32.000000,59.000000i)    (85.000000,100.000000i)    (36.000000,64.000000i)    (65.000000,27.000000i)    
       (51.000000,40.000000i)    (36.000000,64.000000i)    (29.000000,32.000000i)    (26.000000,-11.000000i)    
       (75.000000,-1.000000i)    (65.000000,27.000000i)    (26.000000,-11.000000i)    (6.000000,-48.000000i)    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of C is being reference.	
    trans=Notranspose: both A and B are n by k, return alpha*A*Transpose(A)+beta*C in array C;
    trans=Transpose: both A and B are k by n, return alpha*Transpose(A)*A+beta*C in array C;

Return: NULL
*/

void  zsyrk(
MatrixTriangle uplo, 
MatrixTranspose trans,
int n, //the order of C
int k, 
complex alpha,
const complex a[], //the complex matrix of n by k(trans=Notranspose) or k by n(trans=Transpose)
int tda, //the last dimension of A
complex beta,
complex c[], //Input: the complex symmetric matrix C; Output: matrix C after update
int tdc //the last dimension of C
);

/** f06zwc
       Perform a rank-2k update of a complex symmetric matrix 
Example:
	This example perform a rank-2k update of a complex sysmetric matrix C
   

void test_f06zwc()    //test for function zsyr2k(f06zwc)
{
   int nCols=4;
   int nRows=4;

   matrix<complex>  mat_a={{1,8+3i,2+3i,3+4i},
                           {8-3i,3,5+2i,4+4i},
                           {2-3i,5-2i,2,5+4i},
                           {3-4i,4-4i,5-4i,2}};
                           
   matrix<complex> mat_b={{1+2i,1,1,1},{1,1-4i,1,1},{1,1+1i,1,1},{1-4i,1,1,1}}; 
   
   matrix<complex> mat_c(nRows,nRows);
   mat_c=0;

   zsyr2k(UpperTriangle,NoTranspose,nRows,nCols,1,
                mat_a,nCols,mat_b,nCols,1,mat_c,nRows);


   printf("Perform a rank-2k update of a complex sysmmetric matrix C=C+A*transpose(B)+B*transpose(A):\n");

   int i,j;
   for(i=0;i<nRows;i++)
    	for(j=0;j<i;j++)
    	{
    		mat_c[i][j].m_re=mat_c[j][i].m_re;
    		mat_c[i][j].m_im=mat_c[j][i].m_im;
    	}
   
   printf("A=\n");                     //output general complex matrix mat_a
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_a[i][j].m_re,mat_a[i][j].m_im);
       printf("\n");
   }

   printf("B=\n");                     //output general complex matrix mat_b
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nCols;j++)
          printf("(%f,%fi)    ",mat_b[i][j].m_re,mat_b[i][j].m_im);
       printf("\n");
   }
   
   printf("C=\n");                     //output complex symmetric matrix mat_c
   for(i=0;i<nRows;i++)
   {
       printf("       ");
       for(j=0;j<nRows;j++)
          printf("(%f,%fi)    ",mat_c[i][j].m_re,mat_c[i][j].m_im);
       printf("\n");
   }   
}

   The output is as following:
Perform a rank-2k update of a complex sysmmetric matrix C=C+A*transpose(B)+B*transpose(A):
A=
       (1.000000,0.000000i)    (8.000000,3.000000i)    (2.000000,3.000000i)    (3.000000,4.000000i)    
       (8.000000,-3.000000i)    (3.000000,0.000000i)    (5.000000,2.000000i)    (4.000000,4.000000i)    
       (2.000000,-3.000000i)    (5.000000,-2.000000i)    (2.000000,0.000000i)    (5.000000,4.000000i)    
       (3.000000,-4.000000i)    (4.000000,-4.000000i)    (5.000000,-4.000000i)    (2.000000,0.000000i)    
B=
       (1.000000,2.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,-4.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,0.000000i)    (1.000000,1.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
       (1.000000,-4.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    (1.000000,0.000000i)    
C=
       (28.000000,24.000000i)    (52.000000,-3.000000i)    (31.000000,21.000000i)    (36.000000,0.000000i)    
       (52.000000,-3.000000i)    (40.000000,-18.000000i)    (26.000000,-15.000000i)    (6.000000,-57.000000i)    
       (31.000000,21.000000i)    (26.000000,-15.000000i)    (32.000000,8.000000i)    (20.000000,-17.000000i)    
       (36.000000,0.000000i)    (6.000000,-57.000000i)    (20.000000,-17.000000i)    (-4.000000,-48.000000i)    

Parameters:
	uplo specify whether the upper (uplo=UpperTriangle) or lower (uplo=LowerTriangle) triangle of C is being reference.	
    trans=Notranspose: both A and B are n by k, return alpha*A*transpose(B)+alpha*B*transpose(A)+beta*C in array C;
    trans=Transpose: both A and B are k by n, return alpha*transpose(A)*B+alpha*transpose(B)*A+beta*C in array C;

Return: NULL
*/
void  zsyr2k(
MatrixTriangle uplo, 
MatrixTranspose trans, 
int n, //the order of C
int k, 
complex alpha,
const complex a[], //the general matrix A
int tda, //the last dimension of A
const complex b[], //the general matrix B
int tdb, //the last dimension of B
complex beta, 
complex c[],//Input: the complex symmetric matrix C; Output: matrix C after update
int tdc //the last dimension of C
);
/* end proto */
#endif //!_O_NAG_f01_H