From 2dabc00288b53d6e1bdbf38121a02f141ff73fb5 Mon Sep 17 00:00:00 2001
From: Nicolas Bock <nicolasbock@gmail.com>
Date: Wed, 19 Jul 2017 10:09:27 +0200
Subject: [PATCH 1/2] Remove some unused variables

---
 src/C-interface/dense/bml_adjungate_triangle_dense_typed.c  | 2 --
 src/C-interface/dense/bml_copy_dense_typed.c                | 2 --
 src/C-interface/dense/bml_multiply_dense_typed.c            | 4 ----
 src/C-interface/ellsort/bml_add_ellsort_typed.c             | 2 --
 .../ellsort/bml_adjungate_triangle_ellsort_typed.c          | 6 ++----
 src/C-interface/ellsort/bml_copy_ellsort_typed.c            | 6 ------
 6 files changed, 2 insertions(+), 20 deletions(-)

diff --git a/src/C-interface/dense/bml_adjungate_triangle_dense_typed.c b/src/C-interface/dense/bml_adjungate_triangle_dense_typed.c
index 067dc8ca5..9d009c0c2 100644
--- a/src/C-interface/dense/bml_adjungate_triangle_dense_typed.c
+++ b/src/C-interface/dense/bml_adjungate_triangle_dense_typed.c
@@ -29,9 +29,7 @@ void TYPED_FUNC(
     bml_matrix_dense_t * A,
     char *triangle)
 {
-
     int N = A->N;
-    char mychar;
 
     REAL_T *A_matrix = A->matrix;
 
diff --git a/src/C-interface/dense/bml_copy_dense_typed.c b/src/C-interface/dense/bml_copy_dense_typed.c
index e99f5cdb3..69e299748 100644
--- a/src/C-interface/dense/bml_copy_dense_typed.c
+++ b/src/C-interface/dense/bml_copy_dense_typed.c
@@ -77,13 +77,11 @@ void TYPED_FUNC(
     }
 
     // Reorder elements in each row - just change index
-    REAL_T tmp;
 #pragma omp parallel for
     for (int i = 0; i < N; i++)
     {
         for (int j = 0; j < N; j++)
         {
-            tmp = A_matrix[ROWMAJOR(i, j, N, N)];
             A_matrix[ROWMAJOR(i, j, N, N)] =
                 B_matrix[ROWMAJOR(i, perm[j], N, N)];
         }
diff --git a/src/C-interface/dense/bml_multiply_dense_typed.c b/src/C-interface/dense/bml_multiply_dense_typed.c
index add405295..5bbd2756c 100644
--- a/src/C-interface/dense/bml_multiply_dense_typed.c
+++ b/src/C-interface/dense/bml_multiply_dense_typed.c
@@ -44,12 +44,8 @@ void TYPED_FUNC(
     REAL_T alpha_ = (REAL_T) alpha;
     REAL_T beta_ = (REAL_T) beta;
 
-    int myRank = bml_getMyRank();
-
     C_BLAS(GEMM) ("N", "N", &A->N, &A->N, &A->N, &alpha_, B->matrix,
                   &A->N, A->matrix, &A->N, &beta_, C->matrix, &A->N);
-
-//    mkl_thread_free_buffers();
 }
 
 /** Matrix multiply.
diff --git a/src/C-interface/ellsort/bml_add_ellsort_typed.c b/src/C-interface/ellsort/bml_add_ellsort_typed.c
index 6752243a7..f79d1aca1 100644
--- a/src/C-interface/ellsort/bml_add_ellsort_typed.c
+++ b/src/C-interface/ellsort/bml_add_ellsort_typed.c
@@ -154,8 +154,6 @@ double TYPED_FUNC(
     int *B_nnz = B->nnz;
     int *B_index = B->index;
 
-    int ind, ind2;
-
     REAL_T x[N];
     REAL_T y[N];
     REAL_T *A_value = (REAL_T *) A->value;
diff --git a/src/C-interface/ellsort/bml_adjungate_triangle_ellsort_typed.c b/src/C-interface/ellsort/bml_adjungate_triangle_ellsort_typed.c
index 6133e6455..bd885a0f1 100644
--- a/src/C-interface/ellsort/bml_adjungate_triangle_ellsort_typed.c
+++ b/src/C-interface/ellsort/bml_adjungate_triangle_ellsort_typed.c
@@ -31,8 +31,6 @@ void TYPED_FUNC(
     int A_M = A->M;
     int l;
     int ll;
-    int j;
-    char mychar;
     REAL_T *A_value = (REAL_T *) A->value;
     int *A_index = A->index;
     int *A_nnz = A->nnz;
@@ -49,7 +47,7 @@ void TYPED_FUNC(
 #endif
 
 #pragma omp parallel for default(none) shared(A_N,A_M,A_index,A_nnz,A_value,lock) \
-      private(j,l,ll)
+      private(l,ll)
             //    WARNING: Please, check for race conditions ...
 
             for (int i = 0; i < A_N; i++)       // For every row
@@ -91,7 +89,7 @@ void TYPED_FUNC(
 #endif
 
 #pragma omp parallel for default(none) shared(lock,A_N,A_M,A_index,A_nnz,A_value) \
-      private(j,l,ll)
+      private(l,ll)
             //    WARNING: Please, check for race conditions and parallel performances ...
             for (int i = 0; i < A_N; i++)
             {
diff --git a/src/C-interface/ellsort/bml_copy_ellsort_typed.c b/src/C-interface/ellsort/bml_copy_ellsort_typed.c
index 6e0ea9e6b..e6e5b4d7d 100644
--- a/src/C-interface/ellsort/bml_copy_ellsort_typed.c
+++ b/src/C-interface/ellsort/bml_copy_ellsort_typed.c
@@ -30,11 +30,9 @@ bml_matrix_ellsort_t *TYPED_FUNC(
     int M = A->M;
 
     int *A_index = A->index;
-    int *A_nnz = A->nnz;
     REAL_T *A_value = A->value;
 
     int *B_index = B->index;
-    int *B_nnz = B->nnz;
     REAL_T *B_value = B->value;
 
     //    memcpy(B->index, A->index, sizeof(int) * A->N * A->M);
@@ -47,7 +45,6 @@ bml_matrix_ellsort_t *TYPED_FUNC(
                M * sizeof(int));
         memcpy(&B_value[ROWMAJOR(i, 0, N, M)], &A_value[ROWMAJOR(i, 0, N, M)],
                M * sizeof(REAL_T));
-        //      A_nnz[perm[i]] = B_nnz[i];
     }
     bml_copy_domain(A->domain, B->domain);
     bml_copy_domain(A->domain2, B->domain2);
@@ -70,11 +67,9 @@ void TYPED_FUNC(
     int M = A->M;
 
     int *A_index = A->index;
-    int *A_nnz = A->nnz;
     REAL_T *A_value = A->value;
 
     int *B_index = B->index;
-    int *B_nnz = B->nnz;
     REAL_T *B_value = B->value;
     // memcpy(B->index, A->index, sizeof(int) * A->N * A->M);
     memcpy(B->nnz, A->nnz, sizeof(int) * A->N);
@@ -86,7 +81,6 @@ void TYPED_FUNC(
                M * sizeof(int));
         memcpy(&B_value[ROWMAJOR(i, 0, N, M)], &A_value[ROWMAJOR(i, 0, N, M)],
                M * sizeof(REAL_T));
-        //      A_nnz[perm[i]] = B_nnz[i];
     }
     if (A->distribution_mode == B->distribution_mode)
     {

From b93a366260afdd0ac0b0be4a5970e51173c9b7ed Mon Sep 17 00:00:00 2001
From: Nicolas Bock <nicolasbock@gmail.com>
Date: Wed, 19 Jul 2017 10:09:55 +0200
Subject: [PATCH 2/2] Add complex scale factors

This change change the `bml_scale*` APIs so that the scale factor is
matrix type dependent. This means that a complex matrix type will
require a complex scaling factor.

Fixes: #75
---
 src/C-interface/bml_scale.c                   |  6 +-
 src/C-interface/bml_scale.h                   |  8 +-
 src/C-interface/dense/bml_add_dense_typed.c   |  4 +-
 src/C-interface/dense/bml_scale_dense.c       |  6 +-
 src/C-interface/dense/bml_scale_dense.h       | 32 +++----
 src/C-interface/dense/bml_scale_dense_typed.c | 23 ++---
 .../ellpack/bml_normalize_ellpack_typed.c     |  4 +-
 src/C-interface/ellpack/bml_scale_ellpack.c   |  6 +-
 src/C-interface/ellpack/bml_scale_ellpack.h   | 32 +++----
 .../ellpack/bml_scale_ellpack_typed.c         | 23 ++---
 .../ellsort/bml_normalize_ellsort_typed.c     |  4 +-
 src/C-interface/ellsort/bml_scale_ellsort.c   |  6 +-
 src/C-interface/ellsort/bml_scale_ellsort.h   | 32 +++----
 .../ellsort/bml_scale_ellsort_typed.c         | 23 ++---
 src/Fortran-interface/bml_c_interface_m.F90   |  6 +-
 src/Fortran-interface/bml_scale_m.F90         | 84 +++++++++++++++++--
 src/Fortran-interface/bml_types_m.F90         | 38 ---------
 tests/normalize_matrix_typed.c                |  4 +-
 tests/scale_matrix_typed.c                    |  8 +-
 tests/trace_matrix_typed.c                    | 10 +--
 20 files changed, 182 insertions(+), 177 deletions(-)

diff --git a/src/C-interface/bml_scale.c b/src/C-interface/bml_scale.c
index c045e4806..3068ec862 100644
--- a/src/C-interface/bml_scale.c
+++ b/src/C-interface/bml_scale.c
@@ -17,7 +17,7 @@
  */
 bml_matrix_t *
 bml_scale_new(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_t * A)
 {
     bml_matrix_t *B = NULL;
@@ -50,7 +50,7 @@ bml_scale_new(
  */
 void
 bml_scale(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_t * A,
     bml_matrix_t * B)
 {
@@ -80,7 +80,7 @@ bml_scale(
  */
 void
 bml_scale_inplace(
-    const double scale_factor,
+    const void *scale_factor,
     bml_matrix_t * A)
 {
     switch (bml_get_type(A))
diff --git a/src/C-interface/bml_scale.h b/src/C-interface/bml_scale.h
index 9a9d93860..85c9f646b 100644
--- a/src/C-interface/bml_scale.h
+++ b/src/C-interface/bml_scale.h
@@ -5,19 +5,17 @@
 
 #include "bml_types.h"
 
-// Scales A and returns a new B
 bml_matrix_t *bml_scale_new(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_t * A);
 
-// Scales A and returns in existing B
 void bml_scale(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_t * A,
     bml_matrix_t * B);
 
 void bml_scale_inplace(
-    const double scale_factor,
+    const void *scale_factor,
     bml_matrix_t * A);
 
 #endif
diff --git a/src/C-interface/dense/bml_add_dense_typed.c b/src/C-interface/dense/bml_add_dense_typed.c
index 21f26a75a..5adfb7892 100644
--- a/src/C-interface/dense/bml_add_dense_typed.c
+++ b/src/C-interface/dense/bml_add_dense_typed.c
@@ -146,6 +146,8 @@ void TYPED_FUNC(
     const double alpha,
     const double beta)
 {
-    bml_scale_inplace_dense(alpha, A);
+    REAL_T _alpha = (REAL_T) alpha;
+
+    bml_scale_inplace_dense(&_alpha, A);
     bml_add_identity_dense(A, beta);
 }
diff --git a/src/C-interface/dense/bml_scale_dense.c b/src/C-interface/dense/bml_scale_dense.c
index 68a59b157..240c611ad 100644
--- a/src/C-interface/dense/bml_scale_dense.c
+++ b/src/C-interface/dense/bml_scale_dense.c
@@ -20,7 +20,7 @@
  */
 bml_matrix_dense_t *
 bml_scale_dense_new(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_dense_t * A)
 {
     bml_matrix_dense_t *B = NULL;
@@ -57,7 +57,7 @@ bml_scale_dense_new(
  */
 void
 bml_scale_dense(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_dense_t * A,
     bml_matrix_dense_t * B)
 {
@@ -85,7 +85,7 @@ bml_scale_dense(
 
 void
 bml_scale_inplace_dense(
-    const double scale_factor,
+    const void *scale_factor,
     bml_matrix_dense_t * A)
 {
     switch (A->matrix_precision)
diff --git a/src/C-interface/dense/bml_scale_dense.h b/src/C-interface/dense/bml_scale_dense.h
index 64a46ad67..03cb83c1e 100644
--- a/src/C-interface/dense/bml_scale_dense.h
+++ b/src/C-interface/dense/bml_scale_dense.h
@@ -3,69 +3,71 @@
 
 #include "bml_types_dense.h"
 
+#include <complex.h>
+
 bml_matrix_dense_t *bml_scale_dense_new(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_dense_t * A);
 
 bml_matrix_dense_t *bml_scale_dense_new_single_real(
-    const double scale_factor,
+    const float *scale_factor,
     const bml_matrix_dense_t * A);
 
 bml_matrix_dense_t *bml_scale_dense_new_double_real(
-    const double scale_factor,
+    const double *scale_factor,
     const bml_matrix_dense_t * A);
 
 bml_matrix_dense_t *bml_scale_dense_new_single_complex(
-    const double scale_factor,
+    const float complex * scale_factor,
     const bml_matrix_dense_t * A);
 
 bml_matrix_dense_t *bml_scale_dense_new_double_complex(
-    const double scale_factor,
+    const double complex * scale_factor,
     const bml_matrix_dense_t * A);
 
 void bml_scale_dense(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_dense_t * A,
     bml_matrix_dense_t * B);
 
 void bml_scale_dense_single_real(
-    const double scale_factor,
+    const float *scale_factor,
     const bml_matrix_dense_t * A,
     bml_matrix_dense_t * B);
 
 void bml_scale_dense_double_real(
-    const double scale_factor,
+    const double *scale_factor,
     const bml_matrix_dense_t * A,
     bml_matrix_dense_t * B);
 
 void bml_scale_dense_single_complex(
-    const double scale_factor,
+    const float complex * scale_factor,
     const bml_matrix_dense_t * A,
     bml_matrix_dense_t * B);
 
 void bml_scale_dense_double_complex(
-    const double scale_factor,
+    const double complex * scale_factor,
     const bml_matrix_dense_t * A,
     bml_matrix_dense_t * B);
 
 void bml_scale_inplace_dense(
-    const double scale_factor,
+    const void *scale_factor,
     bml_matrix_dense_t * A);
 
 void bml_scale_inplace_dense_single_real(
-    const double scale_factor,
+    const float *scale_factor,
     bml_matrix_dense_t * A);
 
 void bml_scale_inplace_dense_double_real(
-    const double scale_factor,
+    const double *scale_factor,
     bml_matrix_dense_t * A);
 
 void bml_scale_inplace_dense_single_complex(
-    const double scale_factor,
+    const float complex * scale_factor,
     bml_matrix_dense_t * A);
 
 void bml_scale_inplace_dense_double_complex(
-    const double scale_factor,
+    const double complex * scale_factor,
     bml_matrix_dense_t * A);
 
 #endif
diff --git a/src/C-interface/dense/bml_scale_dense_typed.c b/src/C-interface/dense/bml_scale_dense_typed.c
index f5e9c75ba..980f704d8 100644
--- a/src/C-interface/dense/bml_scale_dense_typed.c
+++ b/src/C-interface/dense/bml_scale_dense_typed.c
@@ -26,21 +26,17 @@
  */
 bml_matrix_dense_t *TYPED_FUNC(
     bml_scale_dense_new) (
-    const double scale_factor,
+    const REAL_T * scale_factor,
     const bml_matrix_dense_t * A)
 {
-    REAL_T sfactor = (REAL_T) scale_factor;
-
     bml_matrix_dense_t *B = TYPED_FUNC(bml_copy_dense_new) (A);
     REAL_T *B_matrix = B->matrix;
     int myRank = bml_getMyRank();
-    //int nElems = B->N * B->N;
     int nElems = B->domain->localRowExtent[myRank] * B->N;
     int startIndex = B->domain->localDispl[myRank];
     int inc = 1;
 
-    C_BLAS(SCAL) (&nElems, &sfactor, &(B_matrix[startIndex]), &inc);
-    //C_BLAS(SCAL) (&nElems, &sfactor, B->matrix, &inc);
+    C_BLAS(SCAL) (&nElems, scale_factor, &(B_matrix[startIndex]), &inc);
 
     return B;
 }
@@ -54,40 +50,33 @@ bml_matrix_dense_t *TYPED_FUNC(
  */
 void TYPED_FUNC(
     bml_scale_dense) (
-    const double scale_factor,
+    const REAL_T * scale_factor,
     const bml_matrix_dense_t * A,
     bml_matrix_dense_t * B)
 {
-    REAL_T sfactor = scale_factor;
-
     if (A != B)
         TYPED_FUNC(bml_copy_dense) (A, B);
 
     REAL_T *B_matrix = B->matrix;
     int myRank = bml_getMyRank();
-    //int nElems = B->N * B->N;
     int nElems = B->domain->localRowExtent[myRank] * B->N;
     int startIndex = B->domain->localDispl[myRank];
     int inc = 1;
 
-    C_BLAS(SCAL) (&nElems, &sfactor, &(B_matrix[startIndex]), &inc);
-    //C_BLAS(SCAL) (&nElems, &sfactor, B->matrix, &inc);
+    C_BLAS(SCAL) (&nElems, scale_factor, &(B_matrix[startIndex]), &inc);
 }
 
 void TYPED_FUNC(
     bml_scale_inplace_dense) (
-    const double scale_factor,
+    const REAL_T * scale_factor,
     bml_matrix_dense_t * A)
 {
-    REAL_T scale_factor_ = (REAL_T) scale_factor;
     REAL_T *A_matrix = A->matrix;
     int myRank = bml_getMyRank();
-    //int number_elements = A->N * A->N;
     int number_elements = A->domain->localRowExtent[myRank] * A->N;
     int startIndex = A->domain->localDispl[myRank];
     int inc = 1;
 
-    C_BLAS(SCAL) (&number_elements, &scale_factor_, &(A_matrix[startIndex]),
+    C_BLAS(SCAL) (&number_elements, scale_factor, &(A_matrix[startIndex]),
                   &inc);
-    //C_BLAS(SCAL) (&number_elements, &scale_factor_, A->matrix, &inc);
 }
diff --git a/src/C-interface/ellpack/bml_normalize_ellpack_typed.c b/src/C-interface/ellpack/bml_normalize_ellpack_typed.c
index 4ed0b2c35..9205777d0 100644
--- a/src/C-interface/ellpack/bml_normalize_ellpack_typed.c
+++ b/src/C-interface/ellpack/bml_normalize_ellpack_typed.c
@@ -36,10 +36,10 @@ void TYPED_FUNC(
 {
     double maxminusmin = maxeval - mineval;
     double gershfact = maxeval / maxminusmin;
-    double scalar = (double) -1.0 / maxminusmin;
+    REAL_T scalar = (REAL_T) - 1.0 / maxminusmin;
     double threshold = 0.0;
 
-    bml_scale_inplace_ellpack(scalar, A);
+    bml_scale_inplace_ellpack(&scalar, A);
     bml_add_identity_ellpack(A, gershfact, threshold);
 }
 
diff --git a/src/C-interface/ellpack/bml_scale_ellpack.c b/src/C-interface/ellpack/bml_scale_ellpack.c
index 683328d9a..8b0b5ab6d 100644
--- a/src/C-interface/ellpack/bml_scale_ellpack.c
+++ b/src/C-interface/ellpack/bml_scale_ellpack.c
@@ -16,7 +16,7 @@
  */
 bml_matrix_ellpack_t *
 bml_scale_ellpack_new(
-    const double scale_factor,
+    const void * scale_factor,
     const bml_matrix_ellpack_t * A)
 {
     bml_matrix_ellpack_t *B = NULL;
@@ -51,7 +51,7 @@ bml_scale_ellpack_new(
  */
 void
 bml_scale_ellpack(
-    const double scale_factor,
+    const void * scale_factor,
     const bml_matrix_ellpack_t * A,
     const bml_matrix_ellpack_t * B)
 {
@@ -77,7 +77,7 @@ bml_scale_ellpack(
 
 void
 bml_scale_inplace_ellpack(
-    const double scale_factor,
+    const void * scale_factor,
     bml_matrix_ellpack_t * A)
 {
     switch (A->matrix_precision)
diff --git a/src/C-interface/ellpack/bml_scale_ellpack.h b/src/C-interface/ellpack/bml_scale_ellpack.h
index 29e15494b..f1fff092b 100644
--- a/src/C-interface/ellpack/bml_scale_ellpack.h
+++ b/src/C-interface/ellpack/bml_scale_ellpack.h
@@ -3,69 +3,71 @@
 
 #include "bml_types_ellpack.h"
 
+#include <complex.h>
+
 bml_matrix_ellpack_t *bml_scale_ellpack_new(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_ellpack_t * A);
 
 bml_matrix_ellpack_t *bml_scale_ellpack_new_single_real(
-    const double scale_factor,
+    const float *scale_factor,
     const bml_matrix_ellpack_t * A);
 
 bml_matrix_ellpack_t *bml_scale_ellpack_new_double_real(
-    const double scale_factor,
+    const double *scale_factor,
     const bml_matrix_ellpack_t * A);
 
 bml_matrix_ellpack_t *bml_scale_ellpack_new_single_complex(
-    const double scale_factor,
+    const float complex * scale_factor,
     const bml_matrix_ellpack_t * A);
 
 bml_matrix_ellpack_t *bml_scale_ellpack_new_double_complex(
-    const double scale_factor,
+    const double complex * scale_factor,
     const bml_matrix_ellpack_t * A);
 
 void bml_scale_ellpack(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_ellpack_t * A,
     const bml_matrix_ellpack_t * B);
 
 void bml_scale_ellpack_single_real(
-    const double scale_factor,
+    const float *scale_factor,
     const bml_matrix_ellpack_t * A,
     const bml_matrix_ellpack_t * B);
 
 void bml_scale_ellpack_double_real(
-    const double scale_factor,
+    const double *scale_factor,
     const bml_matrix_ellpack_t * A,
     const bml_matrix_ellpack_t * B);
 
 void bml_scale_ellpack_single_complex(
-    const double scale_factor,
+    const float complex *scale_factor,
     const bml_matrix_ellpack_t * A,
     const bml_matrix_ellpack_t * B);
 
 void bml_scale_ellpack_double_complex(
-    const double scale_factor,
+    const double complex *scale_factor,
     const bml_matrix_ellpack_t * A,
     const bml_matrix_ellpack_t * B);
 
 void bml_scale_inplace_ellpack(
-    const double scale_factor,
+    const void *scale_factor,
     bml_matrix_ellpack_t * A);
 
 void bml_scale_inplace_ellpack_single_real(
-    const double scale_factor,
+    const float *scale_factor,
     bml_matrix_ellpack_t * A);
 
 void bml_scale_inplace_ellpack_double_real(
-    const double scale_factor,
+    const double *scale_factor,
     bml_matrix_ellpack_t * A);
 
 void bml_scale_inplace_ellpack_single_complex(
-    const double scale_factor,
+    const float complex *scale_factor,
     bml_matrix_ellpack_t * A);
 
 void bml_scale_inplace_ellpack_double_complex(
-    const double scale_factor,
+    const double complex *scale_factor,
     bml_matrix_ellpack_t * A);
 
 #endif
diff --git a/src/C-interface/ellpack/bml_scale_ellpack_typed.c b/src/C-interface/ellpack/bml_scale_ellpack_typed.c
index 972cb7225..3580e4b70 100644
--- a/src/C-interface/ellpack/bml_scale_ellpack_typed.c
+++ b/src/C-interface/ellpack/bml_scale_ellpack_typed.c
@@ -25,22 +25,18 @@
  */
 bml_matrix_ellpack_t *TYPED_FUNC(
     bml_scale_ellpack_new) (
-    const double scale_factor,
+    const REAL_T * scale_factor,
     const bml_matrix_ellpack_t * A)
 {
-    REAL_T sfactor = (REAL_T) scale_factor;
-
     bml_matrix_ellpack_t *B = TYPED_FUNC(bml_copy_ellpack_new) (A);
 
     REAL_T *B_value = B->value;
     int myRank = bml_getMyRank();
-    //int nElems = B->N * B->M;
     int nElems = B->domain->localRowExtent[myRank] * B->M;
     int startIndex = B->domain->localDispl[myRank];
     int inc = 1;
 
-    C_BLAS(SCAL) (&nElems, &sfactor, &(B_value[startIndex]), &inc);
-    //C_BLAS(SCAL) (&nElems, &sfactor, B->value, &inc);
+    C_BLAS(SCAL) (&nElems, scale_factor, &(B_value[startIndex]), &inc);
 
     return B;
 }
@@ -54,41 +50,34 @@ bml_matrix_ellpack_t *TYPED_FUNC(
  */
 void TYPED_FUNC(
     bml_scale_ellpack) (
-    const double scale_factor,
+    const REAL_T * scale_factor,
     const bml_matrix_ellpack_t * A,
     const bml_matrix_ellpack_t * B)
 {
-    REAL_T sfactor = (REAL_T) scale_factor;
-
     if (A != B)
         TYPED_FUNC(bml_copy_ellpack) (A, B);
 
     REAL_T *B_value = B->value;
     int myRank = bml_getMyRank();
-    //int nElems = B->N * B->M;
     int nElems = B->domain->localRowExtent[myRank] * B->M;
     int startIndex = B->domain->localDispl[myRank];
     int inc = 1;
 
-    C_BLAS(SCAL) (&nElems, &sfactor, &(B_value[startIndex]), &inc);
-    //C_BLAS(SCAL) (&nElems, &sfactor, B->value, &inc);
+    C_BLAS(SCAL) (&nElems, scale_factor, &(B_value[startIndex]), &inc);
 }
 
 void TYPED_FUNC(
     bml_scale_inplace_ellpack) (
-    const double scale_factor,
+    const REAL_T * scale_factor,
     bml_matrix_ellpack_t * A)
 {
     REAL_T *A_value = A->value;
-    REAL_T scale_factor_ = (REAL_T) scale_factor;
 
     int myRank = bml_getMyRank();
-    //int number_elements = A->N * A->M;
     int number_elements = A->domain->localRowExtent[myRank] * A->M;
     int startIndex = A->domain->localDispl[myRank];
     int inc = 1;
 
-    C_BLAS(SCAL) (&number_elements, &scale_factor_, &(A_value[startIndex]),
+    C_BLAS(SCAL) (&number_elements, scale_factor, &(A_value[startIndex]),
                   &inc);
-    //C_BLAS(SCAL) (&number_elements, &scale_factor_, A->value, &inc);
 }
diff --git a/src/C-interface/ellsort/bml_normalize_ellsort_typed.c b/src/C-interface/ellsort/bml_normalize_ellsort_typed.c
index dc174ed97..d59f7bc07 100644
--- a/src/C-interface/ellsort/bml_normalize_ellsort_typed.c
+++ b/src/C-interface/ellsort/bml_normalize_ellsort_typed.c
@@ -36,10 +36,10 @@ void TYPED_FUNC(
 {
     double maxminusmin = maxeval - mineval;
     double gershfact = maxeval / maxminusmin;
-    double scalar = (double) -1.0 / maxminusmin;
+    REAL_T scalar = (REAL_T) - 1.0 / maxminusmin;
     double threshold = 0.0;
 
-    bml_scale_inplace_ellsort(scalar, A);
+    bml_scale_inplace_ellsort(&scalar, A);
     bml_add_identity_ellsort(A, gershfact, threshold);
 }
 
diff --git a/src/C-interface/ellsort/bml_scale_ellsort.c b/src/C-interface/ellsort/bml_scale_ellsort.c
index 0dbc7d2ae..3138ac8c6 100644
--- a/src/C-interface/ellsort/bml_scale_ellsort.c
+++ b/src/C-interface/ellsort/bml_scale_ellsort.c
@@ -16,7 +16,7 @@
  */
 bml_matrix_ellsort_t *
 bml_scale_ellsort_new(
-    const double scale_factor,
+    const void * scale_factor,
     const bml_matrix_ellsort_t * A)
 {
     bml_matrix_ellsort_t *B = NULL;
@@ -51,7 +51,7 @@ bml_scale_ellsort_new(
  */
 void
 bml_scale_ellsort(
-    const double scale_factor,
+    const void * scale_factor,
     const bml_matrix_ellsort_t * A,
     const bml_matrix_ellsort_t * B)
 {
@@ -77,7 +77,7 @@ bml_scale_ellsort(
 
 void
 bml_scale_inplace_ellsort(
-    const double scale_factor,
+    const void * scale_factor,
     bml_matrix_ellsort_t * A)
 {
     switch (A->matrix_precision)
diff --git a/src/C-interface/ellsort/bml_scale_ellsort.h b/src/C-interface/ellsort/bml_scale_ellsort.h
index 520c0e8e8..4d717c7f5 100644
--- a/src/C-interface/ellsort/bml_scale_ellsort.h
+++ b/src/C-interface/ellsort/bml_scale_ellsort.h
@@ -3,69 +3,71 @@
 
 #include "bml_types_ellsort.h"
 
+#include <complex.h>
+
 bml_matrix_ellsort_t *bml_scale_ellsort_new(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_ellsort_t * A);
 
 bml_matrix_ellsort_t *bml_scale_ellsort_new_single_real(
-    const double scale_factor,
+    const float * scale_factor,
     const bml_matrix_ellsort_t * A);
 
 bml_matrix_ellsort_t *bml_scale_ellsort_new_double_real(
-    const double scale_factor,
+    const double *scale_factor,
     const bml_matrix_ellsort_t * A);
 
 bml_matrix_ellsort_t *bml_scale_ellsort_new_single_complex(
-    const double scale_factor,
+    const float complex * scale_factor,
     const bml_matrix_ellsort_t * A);
 
 bml_matrix_ellsort_t *bml_scale_ellsort_new_double_complex(
-    const double scale_factor,
+    const double complex *scale_factor,
     const bml_matrix_ellsort_t * A);
 
 void bml_scale_ellsort(
-    const double scale_factor,
+    const void *scale_factor,
     const bml_matrix_ellsort_t * A,
     const bml_matrix_ellsort_t * B);
 
 void bml_scale_ellsort_single_real(
-    const double scale_factor,
+    const float * scale_factor,
     const bml_matrix_ellsort_t * A,
     const bml_matrix_ellsort_t * B);
 
 void bml_scale_ellsort_double_real(
-    const double scale_factor,
+    const double* scale_factor,
     const bml_matrix_ellsort_t * A,
     const bml_matrix_ellsort_t * B);
 
 void bml_scale_ellsort_single_complex(
-    const double scale_factor,
+    const float complex * scale_factor,
     const bml_matrix_ellsort_t * A,
     const bml_matrix_ellsort_t * B);
 
 void bml_scale_ellsort_double_complex(
-    const double scale_factor,
+    const double complex *scale_factor,
     const bml_matrix_ellsort_t * A,
     const bml_matrix_ellsort_t * B);
 
 void bml_scale_inplace_ellsort(
-    const double scale_factor,
+    const void *scale_factor,
     bml_matrix_ellsort_t * A);
 
 void bml_scale_inplace_ellsort_single_real(
-    const double scale_factor,
+    const float * scale_factor,
     bml_matrix_ellsort_t * A);
 
 void bml_scale_inplace_ellsort_double_real(
-    const double scale_factor,
+    const double* scale_factor,
     bml_matrix_ellsort_t * A);
 
 void bml_scale_inplace_ellsort_single_complex(
-    const double scale_factor,
+    const float complex * scale_factor,
     bml_matrix_ellsort_t * A);
 
 void bml_scale_inplace_ellsort_double_complex(
-    const double scale_factor,
+    const double complex *scale_factor,
     bml_matrix_ellsort_t * A);
 
 #endif
diff --git a/src/C-interface/ellsort/bml_scale_ellsort_typed.c b/src/C-interface/ellsort/bml_scale_ellsort_typed.c
index c33d5f475..e042db875 100644
--- a/src/C-interface/ellsort/bml_scale_ellsort_typed.c
+++ b/src/C-interface/ellsort/bml_scale_ellsort_typed.c
@@ -25,22 +25,18 @@
  */
 bml_matrix_ellsort_t *TYPED_FUNC(
     bml_scale_ellsort_new) (
-    const double scale_factor,
+    const REAL_T * scale_factor,
     const bml_matrix_ellsort_t * A)
 {
-    REAL_T sfactor = (REAL_T) scale_factor;
-
     bml_matrix_ellsort_t *B = TYPED_FUNC(bml_copy_ellsort_new) (A);
 
     REAL_T *B_value = B->value;
     int myRank = bml_getMyRank();
-    //int nElems = B->N * B->M;
     int nElems = B->domain->localRowExtent[myRank] * B->M;
     int startIndex = B->domain->localDispl[myRank];
     int inc = 1;
 
-    C_BLAS(SCAL) (&nElems, &sfactor, &(B_value[startIndex]), &inc);
-    //C_BLAS(SCAL) (&nElems, &sfactor, B->value, &inc);
+    C_BLAS(SCAL) (&nElems, scale_factor, &(B_value[startIndex]), &inc);
 
     return B;
 }
@@ -54,41 +50,34 @@ bml_matrix_ellsort_t *TYPED_FUNC(
  */
 void TYPED_FUNC(
     bml_scale_ellsort) (
-    const double scale_factor,
+    const REAL_T * scale_factor,
     const bml_matrix_ellsort_t * A,
     const bml_matrix_ellsort_t * B)
 {
-    REAL_T sfactor = (REAL_T) scale_factor;
-
     if (A != B)
         TYPED_FUNC(bml_copy_ellsort) (A, B);
 
     REAL_T *B_value = B->value;
     int myRank = bml_getMyRank();
-    //int nElems = B->N * B->M;
     int nElems = B->domain->localRowExtent[myRank] * B->M;
     int startIndex = B->domain->localDispl[myRank];
     int inc = 1;
 
-    C_BLAS(SCAL) (&nElems, &sfactor, &(B_value[startIndex]), &inc);
-    //C_BLAS(SCAL) (&nElems, &sfactor, B->value, &inc);
+    C_BLAS(SCAL) (&nElems, scale_factor, &(B_value[startIndex]), &inc);
 }
 
 void TYPED_FUNC(
     bml_scale_inplace_ellsort) (
-    const double scale_factor,
+    const REAL_T * scale_factor,
     bml_matrix_ellsort_t * A)
 {
     REAL_T *A_value = A->value;
-    REAL_T scale_factor_ = (REAL_T) scale_factor;
 
     int myRank = bml_getMyRank();
-    //int number_elements = A->N * A->M;
     int number_elements = A->domain->localRowExtent[myRank] * A->M;
     int startIndex = A->domain->localDispl[myRank];
     int inc = 1;
 
-    C_BLAS(SCAL) (&number_elements, &scale_factor_, &(A_value[startIndex]),
+    C_BLAS(SCAL) (&number_elements, scale_factor, &(A_value[startIndex]),
                   &inc);
-    //C_BLAS(SCAL) (&number_elements, &scale_factor_, A->value, &inc);
 }
diff --git a/src/Fortran-interface/bml_c_interface_m.F90 b/src/Fortran-interface/bml_c_interface_m.F90
index 88ac9ca4f..dceca1b78 100644
--- a/src/Fortran-interface/bml_c_interface_m.F90
+++ b/src/Fortran-interface/bml_c_interface_m.F90
@@ -363,14 +363,14 @@ end function bml_random_matrix_C
 
     subroutine bml_scale_C(alpha, a, b) bind(C, name="bml_scale")
       import :: C_PTR, C_DOUBLE
-      real(C_DOUBLE), value, intent(in) :: alpha
-      type(C_PTR), value :: a
+      type(C_PTR), value, intent(in) :: alpha
+      type(C_PTR), value, intent(in) :: a
       type(C_PTR), value :: b
     end subroutine bml_scale_C
 
     subroutine bml_scale_inplace_C(alpha, a) bind(C, name="bml_scale_inplace")
       import :: C_PTR, C_DOUBLE
-      real(C_DOUBLE), value, intent(in) :: alpha
+      type(C_PTR), value, intent(in) :: alpha
       type(C_PTR), value :: a
     end subroutine bml_scale_inplace_C
 
diff --git a/src/Fortran-interface/bml_scale_m.F90 b/src/Fortran-interface/bml_scale_m.F90
index 8f24907ba..c24f50ab6 100644
--- a/src/Fortran-interface/bml_scale_m.F90
+++ b/src/Fortran-interface/bml_scale_m.F90
@@ -10,7 +10,10 @@ module bml_scale_m
   !> Scale a matrix.
   interface bml_scale
      module procedure scale_one
-     module procedure scale_two
+     module procedure scale_two_single_real
+     module procedure scale_two_double_real
+     module procedure scale_two_single_complex
+     module procedure scale_two_double_complex
   end interface bml_scale
 
   public :: bml_scale
@@ -25,10 +28,10 @@ module bml_scale_m
   !! \param a The matrix
   subroutine scale_one(alpha, a)
 
-    real(C_DOUBLE), intent(in) :: alpha
+    real(C_DOUBLE), target, intent(in) :: alpha
     type(bml_matrix_t), intent(inout) :: a
 
-    call bml_scale_inplace_C(alpha, a%ptr)
+    call bml_scale_inplace_C(c_loc(alpha), a%ptr)
 
   end subroutine scale_one
 
@@ -39,14 +42,81 @@ end subroutine scale_one
   !! \param alpha The factor
   !! \param a The matrix
   !! \param c The matrix
-  subroutine scale_two(alpha, a, c)
+  subroutine scale_two_single_real(alpha, a, c)
 
-    real(C_DOUBLE), intent(in) :: alpha
+    use bml_introspection_m
+
+    real(C_FLOAT), target, intent(in) :: alpha
+    type(bml_matrix_t), intent(in) :: a
+    type(bml_matrix_t), intent(inout) :: c
+
+    integer :: prec_a, prec_c
+
+    call bml_scale_C(c_loc(alpha), a%ptr, c%ptr)
+
+  end subroutine scale_two_single_real
+
+  !> Scale a bml matrix.
+  !!
+  !! \f$ C \leftarrow \alpha A \f$
+  !!
+  !! \param alpha The factor
+  !! \param a The matrix
+  !! \param c The matrix
+  subroutine scale_two_double_real(alpha, a, c)
+
+    use bml_introspection_m
+
+    real(C_DOUBLE), target, intent(in) :: alpha
     type(bml_matrix_t), intent(in) :: a
     type(bml_matrix_t), intent(inout) :: c
 
-    call bml_scale_C(alpha, a%ptr, c%ptr)
+    integer :: prec_a, prec_c
+
+    call bml_scale_C(c_loc(alpha), a%ptr, c%ptr)
+
+  end subroutine scale_two_double_real
+
+  !> Scale a bml matrix.
+  !!
+  !! \f$ C \leftarrow \alpha A \f$
+  !!
+  !! \param alpha The factor
+  !! \param a The matrix
+  !! \param c The matrix
+  subroutine scale_two_single_complex(alpha, a, c)
+
+    use bml_introspection_m
+
+    complex(C_FLOAT_COMPLEX), target, intent(in) :: alpha
+    type(bml_matrix_t), intent(in) :: a
+    type(bml_matrix_t), intent(inout) :: c
+
+    integer :: prec_a, prec_c
+
+    call bml_scale_C(c_loc(alpha), a%ptr, c%ptr)
+
+  end subroutine scale_two_single_complex
+
+  !> Scale a bml matrix.
+  !!
+  !! \f$ C \leftarrow \alpha A \f$
+  !!
+  !! \param alpha The factor
+  !! \param a The matrix
+  !! \param c The matrix
+  subroutine scale_two_double_complex(alpha, a, c)
+
+    use bml_introspection_m
+
+    complex(C_DOUBLE_COMPLEX), target, intent(in) :: alpha
+    type(bml_matrix_t), intent(in) :: a
+    type(bml_matrix_t), intent(inout) :: c
+
+    integer :: prec_a, prec_c
+
+    call bml_scale_C(c_loc(alpha), a%ptr, c%ptr)
 
-  end subroutine scale_two
+  end subroutine scale_two_double_complex
 
 end module bml_scale_m
diff --git a/src/Fortran-interface/bml_types_m.F90 b/src/Fortran-interface/bml_types_m.F90
index edbfea7c0..ddf5d0f28 100644
--- a/src/Fortran-interface/bml_types_m.F90
+++ b/src/Fortran-interface/bml_types_m.F90
@@ -80,42 +80,4 @@ subroutine bml_deallocate(a)
 
   end subroutine bml_deallocate
 
-  !subroutine bml_vector_t_assign(this, other)
-  !  class(bml_vector_t), intent(inout) :: this
-  !  class(bml_vector_t), intent(in) :: other
-  !
-  !  print *, "Direct assignment between vectors not implemented"
-  !  error stop
-  !
-  !end subroutine bml_vector_t_assign
-  !
-  !
-  !subroutine destruct_bml_vector_t(this)
-  !  type(bml_vector_t), intent(inout) :: this
-  !
-  !  print *, "DESTRUCTOR for bml_vector not implemented yet."
-  !  print *, "You possibly leak memory here."
-  !  this%ptr = C_NULL_PTR
-  !
-  !end subroutine destruct_bml_vector_t
-  !
-  !
-  !subroutine bml_matrix_t_assign(this, other)
-  !  class(bml_matrix_t), intent(out) :: this
-  !  class(bml_matrix_t), intent(in) :: other
-  !
-  !  if (c_associated(other%ptr)) then
-  !    this%ptr = bml_copy_new_C(other%ptr)
-  !  end if
-  !
-  !end subroutine bml_matrix_t_assign
-  !
-  !
-  !subroutine destruct_bml_matrix_t(this)
-  !  type(bml_matrix_t), intent(inout) :: this
-  !
-  !    call bml_deallocate(this)
-  !
-  !end subroutine destruct_bml_matrix_t
-
 end module bml_types_m
diff --git a/tests/normalize_matrix_typed.c b/tests/normalize_matrix_typed.c
index 0894ae31e..03d5084c9 100644
--- a/tests/normalize_matrix_typed.c
+++ b/tests/normalize_matrix_typed.c
@@ -25,11 +25,11 @@ int TYPED_FUNC(
     REAL_T *A_dense = NULL;
     REAL_T *B_dense = NULL;
 
-    double scale_factor = 2.5;
+    REAL_T scale_factor = 2.5;
     double threshold = 0.0;
 
     A = bml_identity_matrix(matrix_type, matrix_precision, N, M, sequential);
-    bml_scale_inplace(scale_factor, A);
+    bml_scale_inplace(&scale_factor, A);
     A_gbnd = bml_gershgorin(A);
 
     A_dense = bml_export_to_dense(A, dense_row_major);
diff --git a/tests/scale_matrix_typed.c b/tests/scale_matrix_typed.c
index d25a05402..62d8d9e02 100644
--- a/tests/scale_matrix_typed.c
+++ b/tests/scale_matrix_typed.c
@@ -20,14 +20,14 @@ int TYPED_FUNC(
     REAL_T *B_dense = NULL;
     REAL_T *C_dense = NULL;
 
-    double scale_factor = 2.0;
+    REAL_T scale_factor = 2.0;
 
     //A = bml_random_matrix(matrix_type, matrix_precision, N, M, sequential);
     A = bml_identity_matrix(matrix_type, matrix_precision, N, M, sequential);
-    B = bml_scale_new(scale_factor, A);
+    B = bml_scale_new(&scale_factor, A);
     C = bml_zero_matrix(matrix_type, matrix_precision, N, M, sequential);
-    bml_scale(scale_factor, A, C);
-    bml_scale(scale_factor, A, A);
+    bml_scale(&scale_factor, A, C);
+    bml_scale(&scale_factor, A, A);
 
     A_dense = bml_convert_to_dense(A, dense_row_major);
     B_dense = bml_convert_to_dense(B, dense_row_major);
diff --git a/tests/trace_matrix_typed.c b/tests/trace_matrix_typed.c
index f93538923..e5dfdd029 100644
--- a/tests/trace_matrix_typed.c
+++ b/tests/trace_matrix_typed.c
@@ -28,15 +28,15 @@ int TYPED_FUNC(
     REAL_T *C_dense = NULL;
 
     double traceA, traceB, traceC;
-    double scalar = 0.8;
+    REAL_T scalar = 0.8;
 
     double rel_diff;
 
     LOG_DEBUG("rel. tolerance = %e\n", REL_TOL);
 
     A = bml_identity_matrix(matrix_type, matrix_precision, N, M, sequential);
-    B = bml_scale_new(scalar, A);
-    C = bml_scale_new(scalar, B);
+    B = bml_scale_new(&scalar, A);
+    C = bml_scale_new(&scalar, B);
 
     traceA = bml_trace(A);
     traceB = bml_trace(B);
@@ -106,8 +106,8 @@ int TYPED_FUNC(
     A = bml_import_from_dense(matrix_type, matrix_precision, dense_row_major,
                               N, M, A_dense, 0, sequential);
     bml_free_memory(A_dense);
-    B = bml_scale_new(scalar, A);
-    C = bml_scale_new(scalar, B);
+    B = bml_scale_new(&scalar, A);
+    C = bml_scale_new(&scalar, B);
 
     traceA = bml_trace(A);
     traceB = bml_trace(B);