blas3pp.cc

LAPACK++ (Linear Algebra PACKage in C++) is a software library for numerical linear algebra that sol

void Blas_Mat_Mat_Mult(const LaGenMatComplex &A, 
		       const LaGenMatComplex &B, LaGenMatComplex &C, 
		       bool hermit_A, bool hermit_B, 
		       LaComplex _alpha, LaComplex _beta)
{
   char transa = hermit_A ? 'C' : 'N';
   char transb = hermit_B ? 'C' : 'N';
   // m, n, k have to be determined according to op(A), op(B)!
   integer m = hermit_A ? A.size(1) : A.size(0);
   integer k = hermit_A ? A.size(0) : A.size(1);
   integer n = hermit_B ? B.size(0) : B.size(1);
   integer lda = A.gdim(0), ldb = B.gdim(0), ldc = C.gdim(0);
   doublecomplex alpha(_alpha);
   doublecomplex beta(_beta);

   // Check for correct matrix sizes
   if (alpha.r != 0.0 || alpha.i != 0.0)
   {
      assert(m == C.size(0));
      assert(n == C.size(1));
      assert(k == (hermit_B ? B.size(1) : B.size(0)) );
   }

   F77NAME(zgemm)(&transa, &transb, &m, &n, &k,
		  &alpha, &A(0,0), &lda, &B(0,0), &ldb,
		  &beta, &C(0,0), &ldc);
}

void Blas_Mat_Mat_Mult(const LaGenMatComplex &A, 
		       const LaGenMatComplex &B, LaGenMatComplex &C, 
		       LaComplex _alpha, LaComplex _beta)
{
   Blas_Mat_Mat_Mult(A, B, C, false, false, _alpha, _beta);
}

void Blas_Mat_Trans_Mat_Mult(const LaGenMatComplex &A, 
			     const LaGenMatComplex &B, LaGenMatComplex &C, 
			     LaComplex _alpha, LaComplex _beta)
{
   Blas_Mat_Mat_Mult(A, B, C, true, false, _alpha, _beta);
}

void Blas_Mat_Mat_Trans_Mult(const LaGenMatComplex &A, 
			     const LaGenMatComplex &B, LaGenMatComplex &C, 
			     LaComplex _alpha, LaComplex _beta)
{
   Blas_Mat_Mat_Mult(A, B, C, false, true, _alpha, _beta);
}
#endif // LA_COMPLEX_SUPPORT

// ////////////////////////////////////////////////////////////
// Scaling 

template<class matT, class vecT>
void mat_scale(matT& A, vecT& tmpvec,
	       typename matT::value_type scalar)
{
   int M=A.size(1);
   // max column-sum
   if (M==1)
   {
      // only one column
      tmpvec.ref(A);
      Blas_Scale(scalar, tmpvec);
   }
   else
   {
      for (int k=0; k<M; ++k)
      {
	 tmpvec.ref( A.col(k) );
	 Blas_Scale(scalar, tmpvec);
      }
   }
}

void Blas_Scale(double da, LaGenMatDouble &A)
{
   LaVectorDouble T;
   mat_scale(A, T, da);
}
// void Blas_Scale(float da, LaGenMatFloat &A)
// {
//    LaVectorFloat T;
//    mat_scale(A, T, da);
// }
// void Blas_Scale(int da, LaGenMatInt &A)
// {
//    LaVectorInt T;
//    mat_scale(A, T, da);
// }
// void Blas_Scale(long int da, LaGenMatLongInt &A)
// {
//    LaVectorLongInt T;
//    mat_scale(A, T, da);
// }

#ifdef LA_COMPLEX_SUPPORT
void Blas_Scale(COMPLEX s, LaGenMatComplex &A)
{
   LaVectorComplex T;
   mat_scale(A, T, s);
}
#endif

// ////////////////////////////////////////////////////////////
// Matrix norms

template<class matT, class vecT>
double max_col_sum(const matT& A, vecT& tmpvec)
{
   int M=A.size(1);
   // max column-sum
   if (M==1)
   {
      // only one column
      tmpvec.ref(A);
      return Blas_Norm1( tmpvec );
   }
   else
   {
      LaVectorDouble R(M);
      for (int k=0; k<M; ++k)
      {
	 tmpvec.ref( A( LaIndex(), LaIndex(k) ) );
	 R(k) = Blas_Norm1( tmpvec );
      }
      return Blas_Norm_Inf(R);
   }
}

double Blas_Norm1(const LaGenMatDouble &A)
{
   LaVectorDouble T;
   return max_col_sum(A, T);
}

// ////////////////////////////////////////
// max row norm

template<class matT, class vecT>
double max_row_sum(const matT& A, vecT& tmpvec)
{
   int M=A.size(0);
   // max row-sum
   if (M==1)
   {
      // only one row
      tmpvec.ref(A);
      return Blas_Norm1( tmpvec );
   }
   else
   {
      LaVectorDouble R(M);
      for (int k=0; k<M; ++k)
      {
	 tmpvec.ref( A( LaIndex(k), LaIndex() ) );
	 R(k) = Blas_Norm1( tmpvec );
      }
      return Blas_Norm_Inf(R);
   }
}

double Blas_Norm_Inf(const LaGenMatDouble &A)
{
   LaVectorDouble T;
   return max_row_sum(A, T);
}

// ////////////////////////////////////////
// max frobenius norm

template<class matT, class vecT>
double max_fro_sum(const matT& A, vecT& tmpvec)
{
   int M=A.size(1);
   // calculate this column-wise
   if (M==1)
   {
      // only one column
      tmpvec.ref(A);
      return Blas_Norm2( tmpvec );
   }
   else
   {
      LaVectorDouble R(M);
      for (int k=0; k<M; ++k)
      {
	 tmpvec.ref( A( LaIndex(), LaIndex(k) ) );
	 R(k) = Blas_Norm2( tmpvec );
      }
      return Blas_Norm2(R);
   }
}

double Blas_NormF(const LaGenMatDouble &A)
{
   LaVectorDouble T;
   return max_fro_sum(A, T);
}

/** DEPRECATED, use Blas_Norm_Inf instead. */
double Norm_Inf(const LaGenMatDouble &A) 
{ return Blas_Norm_Inf(A); }


#ifdef LA_COMPLEX_SUPPORT

double Blas_Norm1(const LaGenMatComplex &A)
{
   LaVectorComplex T;
   return max_col_sum(A, T);
}
double Blas_Norm_Inf(const LaGenMatComplex &A)
{
   LaVectorComplex T;
   return max_row_sum(A, T);
}
double Blas_NormF(const LaGenMatComplex &A)
{
   LaVectorComplex T;
   return max_fro_sum(A, T);
}
/** DEPRECATED, use Blas_Norm_Inf instead. */
double Norm_Inf(const LaGenMatComplex &A) 
{ return Blas_Norm_Inf(A); }

#endif // LA_COMPLEX_SUPPORT


// ////////////////////////////////////////////////////////////

#ifdef __GNUC__
#  define LA_STD_ABS std::abs
#else
#  define LA_STD_ABS fabs
#endif

double Norm_Inf(const LaBandMatDouble &A)
{
    // integer kl = A.subdiags(), ku = A.superdiags(); 
    integer N=A.size(1);
    integer M=N;

    // slow version

    LaVectorDouble R(M);
    integer i;
    integer j;
    for (i=0; i<M; i++)
    {
        R(i) = 0.0;
        for (j=0; j<N; j++)
            R(i) += 
		LA_STD_ABS (A(i,j));
    }

    double max = R(0);

    // report back largest row sum
    for (i=1; i<M; i++)
        if (R(i) > max) max=R(i);

    return max;
}


double Norm_Inf(const LaSymmMatDouble &S)
{
    integer N = S.size(0); // square matrix

    // slow version

    LaVectorDouble R(N);
    integer i; 
    integer j;      

    for (i=0; i<N; i++)
    {
        R(i) = 0.0;
        for (j=0; j<N; j++)
            R(i) += 
		LA_STD_ABS (S(i,j));
    }
     
    double max = R(0);

    // report back largest row sum
    for (i=1; i<N; i++)
        if (R(i) > max) max=R(i);

    return max;
}


double Norm_Inf(const LaSpdMatDouble &S)
{
    integer N = S.size(0); //SPD matrices are square

    // slow version

    LaVectorDouble R(N);
    integer i; 
    integer j;      

    for (i=0; i<N; i++)
    {
        R(i) = 0.0;
        for (j=0; j<N; j++)
            R(i) += 
		LA_STD_ABS (S(i,j));
    }
     
    double max = R(0);

    // report back largest row sum
    for (i=1; i<N; i++)
        if (R(i) > max) max=R(i);

    return max;
}


double Norm_Inf(const LaSymmTridiagMatDouble &S)
{
    integer N = S.size();   // S is square
    LaVectorDouble R(N);

    R(0) = LA_STD_ABS(S(0,0)) + LA_STD_ABS(S(0,1));

    for (integer i=1; i<N-1; i++)
    {
        R(i) = LA_STD_ABS(S(i,i-1)) + LA_STD_ABS(S(i,i)) + LA_STD_ABS(S(i,i+1));
    }

    R(N-1) = LA_STD_ABS(S(N-1,N-2)) + LA_STD_ABS(S(N-1,N-1));

    return Norm_Inf(R);
}


double Norm_Inf(const LaTridiagMatDouble &T)
{
    integer N = T.size();   // T is square
    LaVectorDouble R(N);

    R(0) = LA_STD_ABS(T(0,0)) + LA_STD_ABS(T(0,1));

    for (int i=1; i<N-1; i++)
    {
        R(i) = LA_STD_ABS(T(i,i-1)) + LA_STD_ABS(T(i,i)) + LA_STD_ABS(T(i,i+1));
    }

    R(N-1) = LA_STD_ABS(T(N-1,N-2)) + LA_STD_ABS(T(N-1,N-1));

    return Norm_Inf(R);
}

// #undef LA_STD_ABS