* Bug 16020 fixed: cov() was slow 23/20923/27
Stéphane MOTTELET [Fri, 5 Apr 2019 10:31:06 +0000 (12:31 +0200)]
http://bugzilla.scilab.org/16020

Includes reduction of memory bandwidth (in gateway and macro),
update of unit test to include complex input and degenerate cases.
NR test checks speedup w.r.t. legacy code and high
memory bandwith code (https://codereview.scilab.org/#/c/20910).

Change-Id: I87841b6325b67ac64c38461f42b04985af1ebd2a

15 files changed:
scilab/Scilab.sln
scilab/modules/elementary_functions/includes/elem_common.h
scilab/modules/statistics/Makefile.am
scilab/modules/statistics/Makefile.in
scilab/modules/statistics/help/en_US/5_multivariate_stats/cov.xml
scilab/modules/statistics/help/ja_JP/5_multivariate_stats/cov.xml [deleted file]
scilab/modules/statistics/includes/statistics_gw.hxx
scilab/modules/statistics/macros/cov.sci
scilab/modules/statistics/sci_gateway/cpp/sci_percent_cov.cpp [new file with mode: 0644]
scilab/modules/statistics/sci_gateway/cpp/statistics_gw.vcxproj
scilab/modules/statistics/sci_gateway/cpp/statistics_gw.vcxproj.filters
scilab/modules/statistics/sci_gateway/statistics_gateway.xml
scilab/modules/statistics/tests/nonreg_tests/bug_16020.tst [new file with mode: 0644]
scilab/modules/statistics/tests/nonreg_tests/bug_16038.tst [new file with mode: 0644]
scilab/modules/statistics/tests/unit_tests/cov.tst

index 8d83495..0d15e80 100644 (file)
@@ -403,6 +403,9 @@ Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "linear_algebra_gw", "module
 EndProject
 Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "statistics_gw", "modules\statistics\sci_gateway\cpp\statistics_gw.vcxproj", "{EAE1009F-B967-43C4-9408-97A37EFA8678}"
        ProjectSection(ProjectDependencies) = postProject
+               {DBC45B0D-6E0A-4107-B284-5A3B0C5BB50D} = {DBC45B0D-6E0A-4107-B284-5A3B0C5BB50D}
+               {0D3FA25B-8116-44EC-A45E-260789DAA3D9} = {0D3FA25B-8116-44EC-A45E-260789DAA3D9}
+               {5B110267-7C18-437C-B87D-DBA2B50729E9} = {5B110267-7C18-437C-B87D-DBA2B50729E9}
                {18F043DA-1DB5-464F-B67D-CF1C23BE7EA0} = {18F043DA-1DB5-464F-B67D-CF1C23BE7EA0}
        EndProjectSection
 EndProject
index eb2aa16..5c9e4e8 100644 (file)
@@ -61,6 +61,9 @@
 
 extern double C2F(dlamch) (const char *_pszCommand, unsigned long int);
 extern double C2F(logp1) (double *_pdblVal);
+
+// dger   performs the rank 1 operation
+extern int C2F(dger) (int *M, int *N, double* alpha, double* DX, int* incx, double* DY, int* incy,  double *DA, int *lda);
 extern int C2F(dgemm) (char *_pstTransA, char *_pstTransB, int *_piN, int *_piM, int *_piK, double *_pdblAlpha, double *_pdblA, int *_piLdA,
                        double *_pdblB, int *_piLdB, double *_pdblBeta, double *_pdblC, int *_piLdC);
 extern int C2F(zgemm) (char *_pstTransA, char *_pstTransB, int *_piN, int *_piM, int *_piK, double *_pdblAlpha, double *_pdblA, int *_piLdA,
@@ -123,6 +126,8 @@ extern int C2F(ztrsm) (char* side, char* uplo, char* trans, char* diag, int* m,
 // dsyrk does a rank k symmetric update
 extern int C2F(dsyrk) (char* uplo, char* trans, int* n, int* k, double* alpha,
                        double* A, int* lda, double* beta, double* B, int* ldb);
+// dsyr   performs the symmetric rank 1 operation
+extern int C2F(dsyr) (char* uplo, int* n, double* alpha, double *x, int* incx, double* A, int* lda);
 // ztrmm multiply by a triangular matrix
 extern int C2F(ztrmm) (char* side, char* uplo, char* trans, char* diag, int* m, int* n, doublecomplex* alphac,
                        doublecomplex* A, int* lda, doublecomplex* B, int* ldb);
index 87ab4c4..3a9a5fa 100644 (file)
@@ -85,6 +85,9 @@ GATEWAY_C_SOURCES = \
        sci_gateway/c/sci_cdfbet.c \
        sci_gateway/c/sci_cdfnbn.c \
        sci_gateway/c/sci_cdfnor.c
+       
+GATEWAY_CPP_SOURCES = \
+       sci_gateway/cpp/sci_percent_cov.cpp
 
 libscistatistics_la_CPPFLAGS = \
     -I$(srcdir)/includes/ \
@@ -101,6 +104,9 @@ libscistatistics_la_CPPFLAGS = \
     -I$(top_srcdir)/modules/output_stream/includes/ \
     -I$(top_srcdir)/modules/localization/includes/ \
     -I$(top_srcdir)/modules/dynamic_link/includes \
+    -I$(top_srcdir)/modules/elementary_functions/includes \
+    -I$(top_srcdir)/modules/io/includes/ \
+    -I$(top_srcdir)/modules/string/includes/ \
     $(AM_CPPFLAGS)
 
 pkglib_LTLIBRARIES = libscistatistics.la
@@ -109,7 +115,7 @@ noinst_LTLIBRARIES = libscistatistics-algo.la
 
 
 libscistatistics_algo_la_SOURCES = $(STATISTICS_C_SOURCES) $(STATISTICS_FORTRAN_SOURCES)
-libscistatistics_la_SOURCES = $(GATEWAY_C_SOURCES)
+libscistatistics_la_SOURCES = $(GATEWAY_C_SOURCES) $(GATEWAY_CPP_SOURCES)
 libscistatistics_algo_la_CPPFLAGS = $(libscistatistics_la_CPPFLAGS)
 
 # For the code check (splint)
index 4a55b38..e594b5f 100644 (file)
@@ -228,7 +228,9 @@ am__objects_3 = sci_gateway/c/libscistatistics_la-sci_cdfchi.lo \
        sci_gateway/c/libscistatistics_la-sci_cdfbet.lo \
        sci_gateway/c/libscistatistics_la-sci_cdfnbn.lo \
        sci_gateway/c/libscistatistics_la-sci_cdfnor.lo
-am_libscistatistics_la_OBJECTS = $(am__objects_3)
+am__objects_4 =  \
+       sci_gateway/cpp/libscistatistics_la-sci_percent_cov.lo
+am_libscistatistics_la_OBJECTS = $(am__objects_3) $(am__objects_4)
 libscistatistics_la_OBJECTS = $(am_libscistatistics_la_OBJECTS)
 AM_V_P = $(am__v_P_@AM_V@)
 am__v_P_ = $(am__v_P_@AM_DEFAULT_V@)
@@ -257,6 +259,7 @@ am__depfiles_remade =  \
        sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdfnor.Plo \
        sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdfpoi.Plo \
        sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdft.Plo \
+       sci_gateway/cpp/$(DEPDIR)/libscistatistics_la-sci_percent_cov.Plo \
        src/c/$(DEPDIR)/libscistatistics_algo_la-CdfBase.Plo \
        src/c/$(DEPDIR)/libscistatistics_algo_la-ipmpar1.Plo \
        src/c/$(DEPDIR)/libscistatistics_algo_la-sci_string_matrix.Plo
@@ -279,6 +282,24 @@ AM_V_CCLD = $(am__v_CCLD_@AM_V@)
 am__v_CCLD_ = $(am__v_CCLD_@AM_DEFAULT_V@)
 am__v_CCLD_0 = @echo "  CCLD    " $@;
 am__v_CCLD_1 = 
+CXXCOMPILE = $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) \
+       $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS)
+LTCXXCOMPILE = $(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) \
+       $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) \
+       $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) \
+       $(AM_CXXFLAGS) $(CXXFLAGS)
+AM_V_CXX = $(am__v_CXX_@AM_V@)
+am__v_CXX_ = $(am__v_CXX_@AM_DEFAULT_V@)
+am__v_CXX_0 = @echo "  CXX     " $@;
+am__v_CXX_1 = 
+CXXLD = $(CXX)
+CXXLINK = $(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) \
+       $(LIBTOOLFLAGS) --mode=link $(CXXLD) $(AM_CXXFLAGS) \
+       $(CXXFLAGS) $(AM_LDFLAGS) $(LDFLAGS) -o $@
+AM_V_CXXLD = $(am__v_CXXLD_@AM_V@)
+am__v_CXXLD_ = $(am__v_CXXLD_@AM_DEFAULT_V@)
+am__v_CXXLD_0 = @echo "  CXXLD   " $@;
+am__v_CXXLD_1 = 
 F77COMPILE = $(F77) $(AM_FFLAGS) $(FFLAGS)
 LTF77COMPILE = $(LIBTOOL) $(AM_V_lt) --tag=F77 $(AM_LIBTOOLFLAGS) \
        $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS)
@@ -703,6 +724,9 @@ GATEWAY_C_SOURCES = \
        sci_gateway/c/sci_cdfnbn.c \
        sci_gateway/c/sci_cdfnor.c
 
+GATEWAY_CPP_SOURCES = \
+       sci_gateway/cpp/sci_percent_cov.cpp
+
 libscistatistics_la_CPPFLAGS = \
     -I$(srcdir)/includes/ \
     -I$(srcdir)/src/c/ \
@@ -718,12 +742,15 @@ libscistatistics_la_CPPFLAGS = \
     -I$(top_srcdir)/modules/output_stream/includes/ \
     -I$(top_srcdir)/modules/localization/includes/ \
     -I$(top_srcdir)/modules/dynamic_link/includes \
+    -I$(top_srcdir)/modules/elementary_functions/includes \
+    -I$(top_srcdir)/modules/io/includes/ \
+    -I$(top_srcdir)/modules/string/includes/ \
     $(AM_CPPFLAGS)
 
 pkglib_LTLIBRARIES = libscistatistics.la
 noinst_LTLIBRARIES = libscistatistics-algo.la
 libscistatistics_algo_la_SOURCES = $(STATISTICS_C_SOURCES) $(STATISTICS_FORTRAN_SOURCES)
-libscistatistics_la_SOURCES = $(GATEWAY_C_SOURCES)
+libscistatistics_la_SOURCES = $(GATEWAY_C_SOURCES) $(GATEWAY_CPP_SOURCES)
 libscistatistics_algo_la_CPPFLAGS = $(libscistatistics_la_CPPFLAGS)
 
 # For the code check (splint)
@@ -821,7 +848,7 @@ HELP_CHAPTERLANG = en_US fr_FR pt_BR
 all: all-am
 
 .SUFFIXES:
-.SUFFIXES: .sci .bin .c .f .lo .o .obj
+.SUFFIXES: .sci .bin .c .cpp .f .lo .o .obj
 $(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am $(top_srcdir)/Makefile.incl.am $(am__configure_deps)
        @for dep in $?; do \
          case '$(am__configure_deps)' in \
@@ -1082,14 +1109,25 @@ sci_gateway/c/libscistatistics_la-sci_cdfnbn.lo:  \
 sci_gateway/c/libscistatistics_la-sci_cdfnor.lo:  \
        sci_gateway/c/$(am__dirstamp) \
        sci_gateway/c/$(DEPDIR)/$(am__dirstamp)
+sci_gateway/cpp/$(am__dirstamp):
+       @$(MKDIR_P) sci_gateway/cpp
+       @: > sci_gateway/cpp/$(am__dirstamp)
+sci_gateway/cpp/$(DEPDIR)/$(am__dirstamp):
+       @$(MKDIR_P) sci_gateway/cpp/$(DEPDIR)
+       @: > sci_gateway/cpp/$(DEPDIR)/$(am__dirstamp)
+sci_gateway/cpp/libscistatistics_la-sci_percent_cov.lo:  \
+       sci_gateway/cpp/$(am__dirstamp) \
+       sci_gateway/cpp/$(DEPDIR)/$(am__dirstamp)
 
 libscistatistics.la: $(libscistatistics_la_OBJECTS) $(libscistatistics_la_DEPENDENCIES) $(EXTRA_libscistatistics_la_DEPENDENCIES) 
-       $(AM_V_CCLD)$(LINK) -rpath $(pkglibdir) $(libscistatistics_la_OBJECTS) $(libscistatistics_la_LIBADD) $(LIBS)
+       $(AM_V_CXXLD)$(CXXLINK) -rpath $(pkglibdir) $(libscistatistics_la_OBJECTS) $(libscistatistics_la_LIBADD) $(LIBS)
 
 mostlyclean-compile:
        -rm -f *.$(OBJEXT)
        -rm -f sci_gateway/c/*.$(OBJEXT)
        -rm -f sci_gateway/c/*.lo
+       -rm -f sci_gateway/cpp/*.$(OBJEXT)
+       -rm -f sci_gateway/cpp/*.lo
        -rm -f src/c/*.$(OBJEXT)
        -rm -f src/c/*.lo
        -rm -f src/dcdflib/*.$(OBJEXT)
@@ -1109,6 +1147,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdfnor.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdfpoi.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdft.Plo@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscistatistics_la-sci_percent_cov.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscistatistics_algo_la-CdfBase.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscistatistics_algo_la-ipmpar1.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscistatistics_algo_la-sci_string_matrix.Plo@am__quote@ # am--include-marker
@@ -1241,6 +1280,37 @@ sci_gateway/c/libscistatistics_la-sci_cdfnor.lo: sci_gateway/c/sci_cdfnor.c
 @AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@  $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscistatistics_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o sci_gateway/c/libscistatistics_la-sci_cdfnor.lo `test -f 'sci_gateway/c/sci_cdfnor.c' || echo '$(srcdir)/'`sci_gateway/c/sci_cdfnor.c
 
+.cpp.o:
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\
+@am__fastdepCXX_TRUE@  $(CXXCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ $< &&\
+@am__fastdepCXX_TRUE@  $(am__mv) $$depbase.Tpo $$depbase.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXXCOMPILE) -c -o $@ $<
+
+.cpp.obj:
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.obj$$||'`;\
+@am__fastdepCXX_TRUE@  $(CXXCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ `$(CYGPATH_W) '$<'` &&\
+@am__fastdepCXX_TRUE@  $(am__mv) $$depbase.Tpo $$depbase.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXXCOMPILE) -c -o $@ `$(CYGPATH_W) '$<'`
+
+.cpp.lo:
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.lo$$||'`;\
+@am__fastdepCXX_TRUE@  $(LTCXXCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ $< &&\
+@am__fastdepCXX_TRUE@  $(am__mv) $$depbase.Tpo $$depbase.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='$<' object='$@' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LTCXXCOMPILE) -c -o $@ $<
+
+sci_gateway/cpp/libscistatistics_la-sci_percent_cov.lo: sci_gateway/cpp/sci_percent_cov.cpp
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscistatistics_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT sci_gateway/cpp/libscistatistics_la-sci_percent_cov.lo -MD -MP -MF sci_gateway/cpp/$(DEPDIR)/libscistatistics_la-sci_percent_cov.Tpo -c -o sci_gateway/cpp/libscistatistics_la-sci_percent_cov.lo `test -f 'sci_gateway/cpp/sci_percent_cov.cpp' || echo '$(srcdir)/'`sci_gateway/cpp/sci_percent_cov.cpp
+@am__fastdepCXX_TRUE@  $(AM_V_at)$(am__mv) sci_gateway/cpp/$(DEPDIR)/libscistatistics_la-sci_percent_cov.Tpo sci_gateway/cpp/$(DEPDIR)/libscistatistics_la-sci_percent_cov.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='sci_gateway/cpp/sci_percent_cov.cpp' object='sci_gateway/cpp/libscistatistics_la-sci_percent_cov.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscistatistics_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o sci_gateway/cpp/libscistatistics_la-sci_percent_cov.lo `test -f 'sci_gateway/cpp/sci_percent_cov.cpp' || echo '$(srcdir)/'`sci_gateway/cpp/sci_percent_cov.cpp
+
 .f.o:
        $(AM_V_F77)$(F77COMPILE) -c -o $@ $<
 
@@ -1256,6 +1326,7 @@ mostlyclean-libtool:
 clean-libtool:
        -rm -rf .libs _libs
        -rm -rf sci_gateway/c/.libs sci_gateway/c/_libs
+       -rm -rf sci_gateway/cpp/.libs sci_gateway/cpp/_libs
        -rm -rf src/c/.libs src/c/_libs
        -rm -rf src/dcdflib/.libs src/dcdflib/_libs
 install-libscistatistics_la_etcDATA: $(libscistatistics_la_etc_DATA)
@@ -1443,6 +1514,8 @@ distclean-generic:
        -test . = "$(srcdir)" || test -z "$(CONFIG_CLEAN_VPATH_FILES)" || rm -f $(CONFIG_CLEAN_VPATH_FILES)
        -rm -f sci_gateway/c/$(DEPDIR)/$(am__dirstamp)
        -rm -f sci_gateway/c/$(am__dirstamp)
+       -rm -f sci_gateway/cpp/$(DEPDIR)/$(am__dirstamp)
+       -rm -f sci_gateway/cpp/$(am__dirstamp)
        -rm -f src/c/$(DEPDIR)/$(am__dirstamp)
        -rm -f src/c/$(am__dirstamp)
        -rm -f src/dcdflib/$(DEPDIR)/$(am__dirstamp)
@@ -1468,6 +1541,7 @@ distclean: distclean-am
        -rm -f sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdfnor.Plo
        -rm -f sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdfpoi.Plo
        -rm -f sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdft.Plo
+       -rm -f sci_gateway/cpp/$(DEPDIR)/libscistatistics_la-sci_percent_cov.Plo
        -rm -f src/c/$(DEPDIR)/libscistatistics_algo_la-CdfBase.Plo
        -rm -f src/c/$(DEPDIR)/libscistatistics_algo_la-ipmpar1.Plo
        -rm -f src/c/$(DEPDIR)/libscistatistics_algo_la-sci_string_matrix.Plo
@@ -1530,6 +1604,7 @@ maintainer-clean: maintainer-clean-am
        -rm -f sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdfnor.Plo
        -rm -f sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdfpoi.Plo
        -rm -f sci_gateway/c/$(DEPDIR)/libscistatistics_la-sci_cdft.Plo
+       -rm -f sci_gateway/cpp/$(DEPDIR)/libscistatistics_la-sci_percent_cov.Plo
        -rm -f src/c/$(DEPDIR)/libscistatistics_algo_la-CdfBase.Plo
        -rm -f src/c/$(DEPDIR)/libscistatistics_algo_la-ipmpar1.Plo
        -rm -f src/c/$(DEPDIR)/libscistatistics_algo_la-sci_string_matrix.Plo
index cb3e8ed..a2b3966 100644 (file)
@@ -5,6 +5,7 @@
  * Copyright (C) 2012-2013 - Michael Baudin
  * Copyright (C) 2009-2010 - DIGITEO - Michael Baudin
  * Copyright (C) 1993 - 1995 - Anders Holtsberg
+ * Copyright (C) 2019 - Stéphane Mottelet
  *
  * Copyright (C) 2012 - 2016 - Scilab Enterprises
  *
@@ -27,7 +28,7 @@
           xmlns:db="http://docbook.org/ns/docbook">
     <refnamediv>
         <refname>cov</refname>
-        <refpurpose>Covariance matrix</refpurpose>
+        <refpurpose>Sample covariance matrix</refpurpose>
     </refnamediv>
     <refsynopsisdiv>
         <title>Syntax</title>
             <varlistentry>
                 <term>x</term>
                 <listitem>
-                    <para> a nobs-by-1 or nobs-by-nvar matrix of doubles</para>
+                    <para> a nobs-by-1 or nobs-by-n matrix of doubles</para>
                 </listitem>
             </varlistentry>
             <varlistentry>
                 <term>y</term>
                 <listitem>
-                    <para> a nobs-by-1 or nobs-by-nvar matrix of doubles</para>
+                    <para> a nobs-by-1 or nobs-by-m matrix of doubles</para>
                 </listitem>
             </varlistentry>
             <varlistentry>
                 <term>C</term>
                 <listitem>
-                    <para> a square matrix of doubles, the empirical covariance</para>
+                    <para>a square matrix of doubles, the empirical covariance or cross-covariance</para>
                 </listitem>
             </varlistentry>
         </variablelist>
         <title>Description</title>
         <para>
             If x is a nobs-by-1 matrix,
-            then <literal>cov(x)</literal> returns the variance of x,
+            then <literal>cov(x)</literal> returns the sample variance of x,
             normalized by nobs-1.
         </para>
         <para>
-            If x is a nobs-by-nvar matrix,
-            then <literal>cov(x)</literal> returns the nvar-by-nvar covariance matrix of the
+            If x is a nobs-by-n matrix,
+            then <literal>cov(x)</literal> returns the n-by-n sample covariance matrix of the
             columns of x, normalized by nobs-1.
-            Here, each column of x is a variable and
+            Here, each column of x is a variable among (1 ... n) and
             each row of x is an observation.
         </para>
         <para>
             If x and y are two nobs-by-1 matrices,
-            then <literal>cov(x, y)</literal> returns the 2-by-2 covariance matrix of x and
+            then <literal>cov(x, y)</literal> returns the 2-by-2 sample covariance matrix of x and
+            y, normalized by nobs-1, where nobs is the number of observations.
+        </para>
+        <para>
+            If x and y are respectively a nobs-by-n and a nobs-by-m matrix
+            then <literal>cov(x, y)</literal> returns the n-by-m sample cross-covariance matrix of x and
             y, normalized by nobs-1, where nobs is the number of observations.
         </para>
         <para>
             <literal>cov(x, 0)</literal> is the same as <literal>cov(x)</literal> and
             <literal>cov(x, y, 0)</literal> is the same as <literal>cov(x, y)</literal>.
             In this case, if the population is from a normal distribution,
-            then C is the best unbiased estimate of the covariance matrix.
+            then C is the best unbiased estimate of the covariance matrix or cross-covariance matrix.
         </para>
         <para>
             <literal>cov(x, 1)</literal> and <literal>cov(x, y, 1)</literal> normalize by nobs.
             observations about their mean.
         </para>
         <para>
-            The covariance of X and Y is defined by:
+            The covariance of two random vectors X and Y is defined by:
         </para>
         <para>
             <latex>
         <para>
             where E is the expectation.
         </para>
-        <para>
-            This function is compatible with Matlab.
-        </para>
-        <para>
-        </para>
     </refsection>
     <refsection>
         <title>Examples</title>
@@ -145,12 +146,24 @@ S = [
 ];
 C = cov(A)
    ]]></programlisting>
+        <programlisting role="example"><![CDATA[
+x = rand(10,3);
+y = sum(x,2);
+C1 = cov(x,y)
+
+// sample cross-covariance matrix is a submatrix of sample covariance matrix
+// of concatenated [x,y]
+
+C1 = cov([x,y])(1:3,4)
+   ]]></programlisting>
+
     </refsection>
     <refsection>
         <title>Bibliography</title>
         <para>
             <ulink url="http://en.wikipedia.org/wiki/Covariance_matrix">Wikipedia: Covariance matrix</ulink>
-        </para>
+            <ulink url="http://en.wikipedia.org/wiki/Cross-covariance_matrix">Wikipedia: Cross-covariance matrix</ulink>
+      </para>
         <para>
             [3] <ulink url="http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc541.htm">NIST/SEMATECH e-Handbook of Statistical Methods, 6.5.4.1. Mean Vector and Covariance Matrix</ulink>
         </para>
@@ -165,6 +178,10 @@ C = cov(A)
                 <revnumber>5.5.0</revnumber>
                 <revdescription>cov function added, to improve mvvacov (deprecated)</revdescription>
             </revision>
+            <revision>
+                <revnumber>6.1</revnumber>
+                <revdescription>cross-covariance computation added</revdescription>
+            </revision>
         </revhistory>
     </refsection>
 </refentry>
diff --git a/scilab/modules/statistics/help/ja_JP/5_multivariate_stats/cov.xml b/scilab/modules/statistics/help/ja_JP/5_multivariate_stats/cov.xml
deleted file mode 100644 (file)
index 6175782..0000000
+++ /dev/null
@@ -1,169 +0,0 @@
-<?xml version="1.0" encoding="UTF-8"?>
-<!--
- *
- * Copyright (C) 2012-2013 - Michael Baudin
- * Copyright (C) 2009-2010 - DIGITEO - Michael Baudin
- * Copyright (C) 1993 - 1995 - Anders Holtsberg
- *
- * Copyright (C) 2012 - 2016 - Scilab Enterprises
- *
- * This file is hereby licensed under the terms of the GNU GPL v2.0,
- * pursuant to article 5.3.4 of the CeCILL v.2.1.
- * This file was originally licensed under the terms of the CeCILL v2.1,
- * and continues to be available under such terms.
- * For more information, see the COPYING file which you should have received
- * along with this program.
- *
- -->
-<refentry version="5.0-subset Scilab" xml:id="cov" xml:lang="ja"
-          xmlns="http://docbook.org/ns/docbook"
-          xmlns:xlink="http://www.w3.org/1999/xlink"
-          xmlns:svg="http://www.w3.org/2000/svg"
-          xmlns:ns3="http://www.w3.org/1999/xhtml"
-          xmlns:mml="http://www.w3.org/1998/Math/MathML"
-          xmlns:scilab="http://www.scilab.org"
-          xmlns:db="http://docbook.org/ns/docbook">
-    <refnamediv>
-        <refname>cov</refname>
-        <refpurpose>共分散行列</refpurpose>
-    </refnamediv>
-    <refsynopsisdiv>
-        <title>呼び出し手順</title>
-        <synopsis>
-            C = cov(x)
-            C = cov(x, 0)
-            C = cov(x, 1)
-            C = cov(x, y)
-            C = cov(x, y, 0)
-            C = cov(x, y, 1)
-        </synopsis>
-    </refsynopsisdiv>
-    <refsection>
-        <title>パラメータ</title>
-        <variablelist>
-            <varlistentry>
-                <term>x</term>
-                <listitem>
-                    <para>nobs x 1 またはnobs x nvar のdouble行列</para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>y</term>
-                <listitem>
-                    <para>nobs x 1 または nobs x nvar のdouble行列</para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
-                <term>C</term>
-                <listitem>
-                    <para>doublesの正方行列, 経験的共分散</para>
-                </listitem>
-            </varlistentry>
-        </variablelist>
-    </refsection>
-    <refsection>
-        <title>説明</title>
-        <para>
-            x が nobs x 1 行列の場合,
-            <literal>cov(x)</literal> はxの共分散を nobs-1 で正規化して返します.
-        </para>
-        <para>
-            x が nobs x nvar 行列の場合,
-            <literal>cov(x)</literal> は xの列の nvar x nvar 共分散行列を
-            nobs - 1で正規化して返します.
-            ここで,xの各列は変数でxの各行は観測値です.
-        </para>
-        <para>
-            x と y が nobs x 1 の行列の場合,
-            <literal>cov(x, y)</literal> は x および y の2 x 2 共分散行列,
-            nobs-1で正規化したものを返します.
-            ただし,nobsは観測値の数です.
-        </para>
-        <para>
-            <literal>cov(x, 0)</literal> は <literal>cov(x)</literal> と同じ,
-            <literal>cov(x, y, 0)</literal> は <literal>cov(x, y)</literal>と同じです.
-            この場合, 母集団が正規分布の場合,
-            Cは共分散行列のバイアスなしの最良の推定値です.
-        </para>
-        <para>
-            <literal>cov(x, 1)</literal> および
-            <literal>cov(x, y, 1)</literal> は nobs で正規化します.
-            この場合, Cは観測量の平均に関する2次モーメント行列です.
-        </para>
-        <para>
-            XおよびYの共分散は次のように定義されます:
-        </para>
-        <para>
-            <latex>
-                Cov(X, Y) = E( (X-E(X)) (Y-E(Y))^T )
-            </latex>
-        </para>
-        <para>
-            ただし, E は期待値です.
-        </para>
-        <para>
-            この関数は Matlabと互換性があります.
-        </para>
-        <para>
-        </para>
-    </refsection>
-    <refsection>
-        <title>例</title>
-        <programlisting role="example"><![CDATA[
-x = [1; 2];
-y = [3; 4];
-C = cov(x, y)
-expected = [0.5, 0.5; 0.5, 0.5];
-C = cov([x, y])
-   ]]></programlisting>
-        <programlisting role="example"><![CDATA[
-x = [230; 181; 165; 150; 97; 192; 181; 189; 172; 170];
-y = [125; 99; 97; 115; 120; 100; 80; 90; 95; 125];
-expected = [
-1152.4556, -88.911111
--88.911111, 244.26667
-];
-C = cov(x, y)
-C = cov([x, y])
-   ]]></programlisting>
-        <programlisting role="example"><![CDATA[
-// Source [3]
-A = [
-4.0 2.0 0.60
-4.2 2.1 0.59
-3.9 2.0 0.58
-4.3 2.1 0.62
-4.1 2.2 0.63
-];
-S = [
-0.025 0.0075 0.00175
-0.0075 0.007 0.00135
-0.00175 0.00135 0.00043
-];
-C = cov(A)
-   ]]></programlisting>
-    </refsection>
-    <refsection>
-        <title>参考文献</title>
-        <para>
-            <ulink url="http://en.wikipedia.org/wiki/Covariance_matrix">Wikipedia: Covariance matrix</ulink>
-        </para>
-        <para>
-            [3] <ulink url="http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc541.htm">NIST/SEMATECH e-Handbook of Statistical Methods, 6.5.4.1. Mean Vector and Covariance Matrix</ulink>
-        </para>
-        <para>
-            "Introduction to probability and statistics for engineers and scientists", Sheldon Ross
-        </para>
-    </refsection>
-    <refsection>
-        <title>履歴</title>
-        <revhistory>
-            <revision>
-                <revnumber>5.5.0</revnumber>
-                <revdescription>mvvacov(廃止予定)を改善するために,
-                    cov 関数が追加されました
-                </revdescription>
-            </revision>
-        </revhistory>
-    </refsection>
-</refentry>
index 988ff5a..080f18e 100644 (file)
 #ifndef __STATISTICS_GW_HXX__
 #define __STATISTICS_GW_HXX__
 
+#include "cpp_gateway_prototype.hxx"
 extern "C"
 {
 #include "dynlib_statistics_gw.h"
 }
 
+CPP_GATEWAY_PROTOTYPE_EXPORT(sci_percent_cov, EXTERN_STATISTICS_GW);
+
 class StatisticsModule
 {
 private :
index 461af42..1470c46 100644 (file)
@@ -1,8 +1,8 @@
-// Copyright (C) 2012-2013 - Michael Baudin
-// Copyright (C) 2009-2010 - DIGITEO - Michael Baudin
 // Copyright (C) 1993 - 1995 - Anders Holtsberg
-//
+// Copyright (C) 2009-2010 - DIGITEO - Michael Baudin
+// Copyright (C) 2012-2013 - Michael Baudin
 // Copyright (C) 2012 - 2016 - Scilab Enterprises
+// Copyright (C) 2019 - Stéphane MOTTELET
 //
 // This file is hereby licensed under the terms of the GNU GPL v2.0,
 // pursuant to article 5.3.4 of the CeCILL v.2.1.
@@ -11,8 +11,8 @@
 // For more information, see the COPYING file which you should have received
 // along with this program.
 
-function C = cov(varargin)
-    // Covariance matrix
+function C = cov(x,y,nrmlztn)
+    // Sample covariance or cross covariance matrix
     //
     // Syntax
     //   C = cov(x)
@@ -24,40 +24,38 @@ function C = cov(varargin)
     //
     // Parameters
     // x: a nobs-by-1 or nobs-by-nvar matrix of doubles
-    // y: a nobs-by-1 or nobs-by-nvar matrix of doubles
-    // C: a square matrix of doubles, the empirical covariance
+    // y: a nobs-by-1 or nobs-by-mvar matrix of doubles
+    // C: a square matrix of doubles, the sample (cross)covariance 
     //
     // Description
-    // If x is a nobs-by-1 matrix,
-    // then cov(x) returns the variance of x,
-    // normalized by nobs-1.
-    //
-    // If x is a nobs-by-nvar matrix,
-    // then cov(x) returns the nvar-by-nvar covariance matrix of the
-    // columns of x, normalized by nobs-1.
-    // Here, each column of x is a variable and
-    // each row of x is an observation.
-    //
-    // If x and y are two nobs-by-1 matrices,
-    // then cov(x, y) returns the 2-by-2 covariance matrix of x and
-    // y, normalized by nobs-1, where nobs is the number of observations.
-    //
+    // If x is a vector of length nobs then cov(x) returns the sample
+    // variance of x, normalized by nobs-1.
+    //
+    // If x  and y are vectors of length nobs whose components are a sample
+    // of nobs observations of two randown variables X and Y ∈ R
+    // cov(x, y) returns the 2-by-2 sample covariance matrix of x and
+    // y, normalized by nobs-1, where nobs is the number of observations, i.e.
+    // an unbiased estimator of the covariance matrix 
+    // E[(X-E[X])^2]         E[(X-E[X])*(Y-E[Y]]
+    // E[(X-E[X])*(Y-E[Y]]   E[(Y-E[Y])^2]
+    //
+    // If x is a nobs-by-nvar matrix whose rows are a sample of nobs
+    // observations of a random vector X ∈ R^nvar
+    // cov(x) returns the nvar-by-nvar sample covariance matrix of X
+    // normalized by nobs-1, i.e. an unbiased estimator of the covariance
+    // matrix E[(X-E[X])(X-E[X])^T]
+    //
+    // If x is a nobs-by-nvar matrix and y a nobs-by-mvar matrix whose
+    // rows are samples of nobs observations of random vectors X ∈ R^nvar
+    // Y ∈ R^mvar
+    // cov(x) returns the nvar-by-mvar sample cross-covariance matrix of X
+    // and Y normalized by nobs-1, i.e. an unbiased estimator of the 
+    // cross-covariance matrix E[(X-E[X])(Y-E[Y])^T]
+    // 
     // cov(x, 0) is the same as cov(x) and
     // cov(x, y, 0) is the same as cov(x, y).
-    // In this case, if the population is from a normal distribution,
-    // then C is the best unbiased estimate of the covariance matrix.
     //
     // cov(x, 1) and cov(x, y, 1) normalize by nobs.
-    // In this case, C is the second moment matrix of the
-    // observations about their mean.
-    //
-    // The covariance of X and Y is defined by:
-    //
-    // Cov(X, Y) = E( (X-E(X)) (Y-E(Y))^T )
-    //
-    // where E is the expectation.
-    //
-    // This function is compatible with Matlab.
     //
     // Examples
     // x = [1; 2];
@@ -98,107 +96,79 @@ function C = cov(varargin)
 
     [lhs, rhs]=argn()
     //
-    if (rhs == 1) then
-        x = varargin(1)
-        //
-        // Check type
-        if (typeof(x) <> "constant")
-            error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
-        end
-        nobs = size(x, "r")
-        r = 1/(nobs-1)
-        A = x
-    elseif (rhs == 2) then
-        //
-        x = varargin(1)
-        y = varargin(2)
+    if rhs ==0
+        error(msprintf(gettext("%s: Wrong number of input argument(s): %d, %d or %d expected.\n"),"cov", 1, 2, 3));
+    elseif rhs == 1
+        y = 0
+    end
+    //
+    // Check type
+    if (typeof(x) <> "constant")
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
+    end
+
+    if (rhs <= 2) then
         //
-        // Check type
-        if (typeof(x) <> "constant")
-            error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
-        end
         if (typeof(y) <> "constant")
             error(msprintf(gettext("%s: Wrong type for input argument #%d: an integer or a real matrix expected.\n"),"cov", 2));
         end
         //
         // Check size
-        nobs = size(x, "r")
-        if (size(y, "*") == 1) then
-            if (y <> 0 & y <> 1)
-                error(msprintf(gettext("%s: Wrong value for input argument #%d: %d or %d expected.\n"),"cov", 2, 0, 1));
+        if isscalar(x) & isscalar(y)
+            if y == 0 || y == 1
+                C = 0;
+            else
+                C = zeros(2,2);
             end
-            if (y == 1) then
-                r = 1/nobs
-                A = x
-            elseif (y == 0) then
-                r = 1/(nobs-1)
-                A = x
-            end
-        else
-            if (size(x) <> [nobs 1]) then
-                error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 1, nobs, 1));
+            return
+        elseif isscalar(y)
+            if y <> 0 & y <> 1
+                error(msprintf(gettext("%s: Wrong value for input argument #%d: %d or %d expected.\n"),"cov", 2, 0, 1));
             end
-            if (size(y) <> [nobs 1]) then
-                error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 2, nobs, 1));
+            //
+            if isvector(x)
+                nobs = length(x)
+                r = 1 / ( nobs - 1 + y)
+                C = r * (norm(x)^2 - nobs*mean(x)^2);
+            else
+                nobs = size(x, "r")
+                r = 1 / (nobs - 1 + y)
+                C = %_cov(x,r);
             end
-            r = 1/(nobs-1)
-            A = [x, y]
-        end
-    elseif (rhs == 3) then
-        //
-        x = varargin(1)
-        y = varargin(2)
-        nrmlztn = varargin(3)
-        //
-        // Check type
-        if (typeof(x) <> "constant")
-            error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1));
+            return
         end
+        nrmlztn = 0;
+    end
+
+    if (rhs == 3)
         if (typeof(y) <> "constant")
             error(msprintf(gettext("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 2));
         end
         if (typeof(nrmlztn) <> "constant")
             error(msprintf(gettext("%s: Wrong type for input argument #%d: an integer expected.\n"),"cov", 3));
         end
-        //
-        // Check size
-        nobs = size(x, "r")
-        if (size(x) <> [nobs 1]) then
-            error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 1, nobs, 1));
-        end
-        if (size(y) <> [nobs 1]) then
-            error(msprintf(gettext("%s: Wrong size for input argument #%d: %dx%d expected.\n"),"cov", 2, nobs, 1));
-        end
-        if (size(nrmlztn, "*") <> 1) then
-            error(msprintf(gettext("%s: Wrong type for input argument #%d: an integer expected.\n"),"cov", 3));
-        end
-        //
-        // Check content
-        if (nrmlztn <> 0 & nrmlztn <> 1)
+        if ( ~isscalar(nrmlztn) || (nrmlztn <> 0 & nrmlztn <> 1))
             error(msprintf(gettext("%s: Wrong value for input argument #%d: %d or %d expected.\n"),"cov", 3, 0, 1));
         end
-        A = [x, y]
-        if (nrmlztn == 1) then
-            r = 1/nobs
-        else
-            r = 1/(nobs-1)
-        end
-    else
-        error(msprintf(gettext("%s: Wrong number of input argument(s): %d, %d or %d expected.\n"),"cov", 1, 2, 3));
     end
-    //
-    // Compute with A in the general case
-    nvar = size(A, "c")
-    nobs = size(A, "r")
-    for i = 1:nvar
-        A(:,i) = A(:,i) - mean(A(:,i))
-    end
-    C = zeros(nvar, nvar)
-    for i = 1:nvar
-        C(i,i) = A(:,i)'*A(:,i)*r
-        for j = i+1:nvar
-            C(i,j) = A(:,i)'*A(:,j)*r
-            C(j,i) = C(i,j)
+
+    if isvector(x) & isvector(y)
+        nobs = length(x)
+        if nobs <> length(y)
+            error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same sizes expected"),"cov", 1, 2));
         end
+        r = 1 / (nobs - 1 + nrmlztn);
+        mx = mean(x);
+        my = mean(y);
+        C = r*[norm(x)^2 - nobs*mx^2, x'*y - nobs*mx*my
+                                   0, norm(y)^2 - nobs*my^2];
+        C(2,1)=C(1,2);
+    elseif size(x,"r") <> size(y,"r")
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same number of rows expected.\n"),"cov", 1, 2));
+    else
+        nobs = size(x,"r")
+        r = 1 / (nobs - 1 + nrmlztn);
+        C = %_cov(x,y,r);
     end
+
 endfunction
diff --git a/scilab/modules/statistics/sci_gateway/cpp/sci_percent_cov.cpp b/scilab/modules/statistics/sci_gateway/cpp/sci_percent_cov.cpp
new file mode 100644 (file)
index 0000000..9f23190
--- /dev/null
@@ -0,0 +1,191 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2019 - Stéphane MOTTELET
+ *
+ * This file is hereby licensed under the terms of the GNU GPL v2.0,
+ * For more information, see the COPYING file which you should have received
+ * along with this program.
+ *
+ */
+
+#include "double.hxx"
+#include "function.hxx"
+#include "statistics_gw.hxx"
+
+extern "C"
+{
+#include "Scierror.h"
+#include "localization.h"
+#include "elem_common.h"
+}
+
+/* ==================================================================== */
+types::Function::ReturnValue sci_percent_cov(types::typed_list& in, int _iRetCount, types::typed_list& out)
+{
+    types::Double* pDblOut = nullptr;
+
+    if (in.size() != 2 && in.size() != 3)
+    {
+        Scierror(77, _("%s: Wrong number of input argument(s): %d or %d expected.\n"), "cov", 2, 3);
+        return types::Function::Error;
+    }
+
+    if (_iRetCount != 1)
+    {
+        Scierror(78, _("%s: Wrong number of output argument(s): %d expected."), "cov", 1);
+        return types::Function::Error;
+    }
+
+    for (int i = 0; i < in.size(); i++)
+    {
+        if (in[i]->isDouble() == false)
+        {
+            Scierror(999, _("%s: Wrong type for input argument #%d: %s expected.\n"), "cov", i + 1, "double");
+            return types::Function::Error;
+        }
+    }
+
+    if (in[in.size() - 1]->getAs<types::Double>()->getSize() != 1)
+    {
+        Scierror(999, _("%s: Wrong dimensions for input argument #%d: A scalar expected.\n"), "cov", in.size());
+        return types::Function::Error;
+    }
+
+    for (int i = 0; i < in.size() - 1; i++)
+    {
+        if (in[i]->getAs<types::Double>()->getDims() > 2)
+        {
+            Scierror(999, _("%s: Wrong dimensions for input argument #%d: A column vector or matrix expected.\n"), "cov", i + 1);
+            return types::Function::Error;
+        }
+    }
+
+    if (in[0]->getAs<types::Double>()->isEmpty() || (in.size() == 3 && in[1]->getAs<types::Double>()->isEmpty()))
+    {
+        pDblOut = types::Double::Empty();
+        out.push_back(pDblOut);
+        return types::Function::OK;
+    }
+
+    double dblr1 = in[in.size() - 1]->getAs<types::Double>()->get(0);
+
+    types::Double* pDblX = in[0]->getAs<types::Double>();
+    int* piXDims = pDblX->getDimsArray();
+    int iNrow = piXDims[0];
+    int iNcolX = piXDims[1];
+    int iOne = 1;
+    double dblOne = 1.0;
+    double dblMinusOne = -1.0;
+    double dblZero = 0.0;
+    double dblr2 = -dblr1 / (double)iNrow;
+
+    // Real and imaginary part of input matrix X1 + %i * X2
+    double* pdblX[3] = { pDblX->get(), pDblX->getImg(), NULL };
+    double* pdblOut[2];
+    // Vector of ones for row-wise sums
+    types::Double* pDblU = new types::Double(iNrow, 1, false);
+    pDblU->setOnes();
+
+    if (in.size() == 3) // covi(X,Y,r) when X and Y are matrices
+    {
+        double dblMinusr1 = -dblr1;
+        double dblMinusr2 = -dblr2;
+        types::Double* pDblY = in[1]->getAs<types::Double>();
+        int* piYDims = pDblY->getDimsArray();
+        int iNcolY = piYDims[1];
+        // Real and imaginary part of input matrix Y = Y1 + %i * Y2
+        double* pdblY[3] = { pDblY->get(), pDblY->getImg(), NULL };
+
+        if (piXDims[0] != piYDims[0])
+        {
+            Scierror(999, _("%s: Incompatible input arguments #%d and #%d: Same number of rows expected.\n"), "cov", 1, 2);
+            return types::Function::Error;
+        }
+
+        pDblOut = new types::Double(iNcolX, iNcolY, pDblX->isComplex() || pDblY->isComplex());
+        pdblOut[0] = pDblOut->get();
+        pdblOut[1] = pDblOut->getImg();
+        // Vector for row-wise sums
+        double *pdblSumX = new double[iNcolX];
+        double *pdblSumY[2] = { new double[iNcolY], new double[iNcolY] };
+
+        pDblOut->setZeros();
+
+        for (int j = 0; pdblY[j] != NULL; j++)
+        {
+            // pre-computation of sum(Yj,'r')
+            C2F(dgemv)((char*)"T", &iNrow, &iNcolY, &dblOne, pdblY[j], &iNrow, pDblU->get(), &iOne, &dblZero, &pdblSumY[j][0], &iOne);
+        }
+        for (int i = 0; pdblX[i] != NULL; i++)
+        {
+            // Computation of sum(Xi,'r')
+            C2F(dgemv)((char*)"T", &iNrow, &iNcolX, &dblOne, pdblX[i], &iNrow, pDblU->get(), &iOne, &dblZero, pdblSumX, &iOne);
+            for (int j = 0; pdblY[j] != NULL; j++)
+            {
+                // +/- r1 * Xi'*Yj
+                C2F(dgemm)((char*)"T", (char*)"N", &iNcolX, &iNcolY, &iNrow, i == 1 && j == 0 ? &dblMinusr1 : &dblr1,
+                    pdblX[i], &iNrow, pdblY[j], &iNrow, &dblOne, pdblOut[i^j], &iNcolX);
+                // +/- r2 * (Xi'*U)*(Yj'*U)'
+                C2F(dger)(&iNcolX, &iNcolY, i == 1 && j == 0 ? &dblMinusr2 : &dblr2,
+                    pdblSumX, &iOne, &pdblSumY[j][0], &iOne, pdblOut[i^j], &iNcolX);
+            }
+        }
+        delete[] pdblSumX;
+        delete[] pdblSumY[0];
+        delete[] pdblSumY[1];
+    }
+    else // (in.size() == 2) covi(X,r)
+    {
+        double* pdblSrc = NULL;
+
+        // Output matrix Re(cov) + %i*Im(cov)
+        pDblOut = new types::Double(iNcolX, iNcolX, pDblX->isComplex());
+        pdblOut[0] = pDblOut->get();
+        pdblOut[1] = pDblOut->getImg();
+
+        // Vector for row-wise sums
+        double* pdblSumX[2] = { new double[iNcolX], new double[iNcolX] };
+
+        pDblOut->setZeros();
+
+        for (int i = 0; pdblX[i] != NULL; i++)
+        {
+            // Re(cov) += r1 * X'*X + r2 * (X'*U)*(U'*X) with X=X1 then X2 if complex input
+            // Computation of sum(Xi,'r')
+            C2F(dgemv)((char*)"T", &iNrow, &iNcolX, &dblOne, pdblX[i], &iNrow, pDblU->get(), &iOne, &dblZero, &pdblSumX[i][0], &iOne);
+            // Below computations are done on upper triangle only
+            C2F(dsyrk)((char*)"U", (char*)"T", &iNcolX, &iNrow, &dblr1, pdblX[i], &iNrow, &dblOne, pdblOut[0], &iNcolX);
+            C2F(dsyr)((char*)"U", &iNcolX, &dblr2, &pdblSumX[i][0], &iOne, pdblOut[0], &iNcolX);
+        }
+
+        if (pDblX->isComplex())
+        {
+            // Im(cov) = r1 * X1'*X2 + r2 * (X1'*U)*(U'*X2) - transpose of previous term
+            // BLAS has no skew symmetric optimized functions, we compute whole first term
+            C2F(dgemm)((char*)"T", (char*)"N", &iNcolX, &iNcolX, &iNrow, &dblr1, pdblX[0], &iNrow, pdblX[1], &iNrow, &dblZero, pdblOut[1], &iNcolX);
+            C2F(dger)(&iNcolX, &iNcolX, &dblr2, &pdblSumX[0][0], &iOne, &pdblSumX[1][0], &iOne, pdblOut[1], &iNcolX);
+            // then substract transpose
+            pdblSrc = pdblOut[1];
+            for (int iLen = iNcolX; iLen > 0; iLen--, pdblSrc += (iNcolX + 1))
+            {
+                C2F(daxpy)(&iLen, &dblMinusOne, pdblSrc, &iOne, pdblSrc, &iNcolX);
+                C2F(dcopy)(&iLen, pdblSrc, &iNcolX, pdblSrc, &iOne);
+                C2F(dscal)(&iLen, &dblMinusOne, pdblSrc, &iOne);
+            }
+        }
+
+        // Real part: copy upper triangle into lower triangle (no BLAS function does this)
+        pdblSrc = pdblOut[0];
+        for (int iLen = iNcolX; iLen > 1; iLen--, pdblSrc += (iNcolX + 1))
+        {
+            C2F(dcopy)(&iLen, pdblSrc, &iNcolX, pdblSrc, &iOne);
+        }
+
+        delete[] pdblSumX[0];
+        delete[] pdblSumX[1];
+    }
+
+    pDblU->killMe();
+    out.push_back(pDblOut);
+    return types::Function::OK;
+}
index 260aad0..e319048 100644 (file)
@@ -88,7 +88,7 @@
       </PrecompiledHeader>
       <WarningLevel>Level3</WarningLevel>
       <DebugInformationFormat>EditAndContinue</DebugInformationFormat>
-      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../core/includes;../../../localization/includes;../../../api_scilab/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../core/includes;../../../localization/includes;../../../api_scilab/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;../../../string/includes;../../../output_stream/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       </PrecompiledHeader>
       <WarningLevel>Level3</WarningLevel>
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
-      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../core/includes;../../../localization/includes;../../../api_scilab/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../core/includes;../../../localization/includes;../../../api_scilab/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;../../../string/includes;../../../output_stream/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <WarningLevel>Level3</WarningLevel>
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
       <MultiProcessorCompilation>true</MultiProcessorCompilation>
-      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../core/includes;../../../localization/includes;../../../api_scilab/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../core/includes;../../../localization/includes;../../../api_scilab/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;../../../string/includes;../../../output_stream/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
       <WarningLevel>Level3</WarningLevel>
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
       <MultiProcessorCompilation>true</MultiProcessorCompilation>
-      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../core/includes;../../../localization/includes;../../../api_scilab/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../core/includes;../../../localization/includes;../../../api_scilab/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;../../../string/includes;../../../output_stream/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <OutputFile>$(SolutionDir)bin\$(ProjectName).dll</OutputFile>
     <ClCompile Include="..\c\sci_cdfnor.c" />
     <ClCompile Include="..\c\sci_cdfpoi.c" />
     <ClCompile Include="..\c\sci_cdft.c" />
+    <ClCompile Include="sci_percent_cov.cpp" />
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="..\..\includes\dynlib_statistics_gw.h" />
     <ClInclude Include="..\..\includes\statistics_gw.hxx" />
   </ItemGroup>
   <ItemGroup>
+    <ProjectReference Include="..\..\..\ast\ast.vcxproj">
+      <Project>{0d3fa25b-8116-44ec-a45e-260789daa3d9}</Project>
+    </ProjectReference>
+    <ProjectReference Include="..\..\..\output_stream\src\c\output_stream.vcxproj">
+      <Project>{a5911cd7-f8e8-440c-a23e-4843a0636f3a}</Project>
+    </ProjectReference>
     <ProjectReference Include="..\..\src\c\statistics.vcxproj">
       <Project>{fe9eb721-b3c1-41d8-b585-3fb3a0083cec}</Project>
     </ProjectReference>
   </ItemGroup>
   <ItemGroup>
+    <Library Include="..\..\..\..\bin\blasplus.lib" />
     <Library Include="..\..\..\..\bin\dcd_f.lib" />
+    <Library Include="..\..\..\..\bin\libintl.lib" />
   </ItemGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
   <ImportGroup Label="ExtensionTargets">
index e6ea9a0..bbbc627 100644 (file)
@@ -48,6 +48,9 @@
     <ClCompile Include="..\c\sci_cdft.c">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="sci_percent_cov.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="..\..\includes\statistics_gw.hxx">
@@ -62,5 +65,7 @@
   </ItemGroup>
   <ItemGroup>
     <Library Include="..\..\..\..\bin\dcd_f.lib" />
+    <Library Include="..\..\..\..\bin\blasplus.lib" />
+    <Library Include="..\..\..\..\bin\libintl.lib" />
   </ItemGroup>
 </Project>
\ No newline at end of file
index 44a8918..fe0aa8d 100644 (file)
@@ -33,4 +33,5 @@
     <gateway name="sci_cdfnor"  function="cdfnor"   type="0" />
     <gateway name="sci_cdfpoi"  function="cdfpoi"   type="0" />
     <gateway name="sci_cdft"    function="cdft"     type="0" />
+    <gateway name="sci_percent_cov"    function="%_cov"     type="1" />
 </module>
diff --git a/scilab/modules/statistics/tests/nonreg_tests/bug_16020.tst b/scilab/modules/statistics/tests/nonreg_tests/bug_16020.tst
new file mode 100644 (file)
index 0000000..67ad223
--- /dev/null
@@ -0,0 +1,49 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2019 - Stéphane MOTTELET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+// <-- Non-regression test for bug 16020 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=16020
+//
+// <-- Short Description -->
+// cov() was slow
+// =============================================================================
+// legacy code
+n=20;
+x=rand(1e6,n);
+for i = 1:n
+    x(:,i) = x(:,i) - mean(x(:,i));
+end
+tic
+c = zeros(n, n);
+for i = 1:n
+    c(i,i) = x(:,i)'*x(:,i);
+    for j = i+1:n
+        c(i,j) = x(:,i)'*x(:,j);
+        c(j,i) = c(i,j);
+    end
+end
+t0 = toc(); // typically t0 = 5 seconds
+
+// large memory bandwidth implementation
+x=rand(1e6,n);
+tic;
+x = center(x,"r");
+c = x'*x;
+t1 = toc(); // typically t0 = 1 seconds
+
+// new implementation (c++ gateway)
+x=rand(1e6,n);
+tic
+c = cov(x);
+t2=toc(); // typically t2 = 0.06 seconds
+
+assert_checktrue(t1/t2 > 4);
+assert_checktrue(t0/t2 > 30);
\ No newline at end of file
diff --git a/scilab/modules/statistics/tests/nonreg_tests/bug_16038.tst b/scilab/modules/statistics/tests/nonreg_tests/bug_16038.tst
new file mode 100644 (file)
index 0000000..043fe52
--- /dev/null
@@ -0,0 +1,26 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2019 - Stéphane MOTTELET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- CLI SHELL MODE -->
+// <-- NO CHECK REF -->
+// <-- Non-regression test for bug 16038 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=16038
+//
+// <-- Short Description -->
+// cov(x,y) size checking of arguments is incorrect
+// =============================================================================
+//
+x = rand(3,1);
+y = rand(4,1);
+message = msprintf(_("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"),"cov",1,2);
+assert_checkerror("cov(x,y)",message);
+x = rand(3,3);
+y = rand(4,4);
+message = msprintf(_("%s: Incompatible input arguments #%d and #%d: Same number of rows expected.\n"),"cov",1,2);
+assert_checkerror("cov(x,y)",message);
index ad8e35b..fc90806 100644 (file)
@@ -1,5 +1,6 @@
 // Copyright (C) 2012-2013 - Michael Baudin
 // Copyright (C) 2010 - INRIA - Michael Baudin
+// Copyright (C) 2019 - Stéphane MOTTELET
 //
 // This file must be used under the terms of the GNU Lesser General Public License license :
 // http://www.gnu.org/copyleft/lesser.html
@@ -7,92 +8,92 @@
 // Run with test_run('statistics', 'cov', ['no_check_error_output']);
 
 // <-- CLI SHELL MODE -->
-// <-- ENGLISH IMPOSED -->
+// <-- NO CHECK REF -->
 
 // Check error
-assert_checkfalse(execstr("cov()"   ,"errcatch") == 0);
-refMsg = msprintf(_("%s: Wrong number of input argument(s): %d, %d or %d expected.\n"),"cov", 1, 2, 3);
+assert_checkfalse(execstr("cov()", "errcatch") =  = 0);
+refMsg = msprintf(_("%s: Wrong number of input argument(s): %d, %d or %d expected.\n"), "cov", 1, 2, 3);
 assert_checkerror("cov()", refMsg);
 
-assert_checkfalse(execstr("cov(""r"")"   ,"errcatch") == 0);
-refMsg = msprintf(_("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1);
+assert_checkfalse(execstr("cov(""r"")", "errcatch") =  = 0);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: a real matrix expected.\n"), "cov", 1);
 assert_checkerror("cov(""r"")", refMsg);
 
-assert_checkfalse(execstr("cov([1;2], ""r"")"   ,"errcatch") == 0);
-refMsg = msprintf(_("%s: Wrong type for input argument #%d: an integer or a real matrix expected.\n"),"cov", 2);
+assert_checkfalse(execstr("cov([1;2], ""r"")", "errcatch") =  = 0);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: an integer or a real matrix expected.\n"), "cov", 2);
 assert_checkerror("cov([1;2], ""r"")", refMsg);
 
-assert_checkfalse(execstr("cov(""r"", [1;2])"   ,"errcatch") == 0);
-refMsg = msprintf(_("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1);
+assert_checkfalse(execstr("cov(""r"", [1;2])", "errcatch") =  = 0);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: a real matrix expected.\n"), "cov", 1);
 assert_checkerror("cov(""r"", [1;2])", refMsg);
 
-assert_checkfalse(execstr("cov(""r"", [1;2], 1)"   ,"errcatch") == 0);
-refMsg = msprintf(_("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 1);
+assert_checkfalse(execstr("cov(""r"", [1;2], 1)", "errcatch") =  = 0);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: a real matrix expected.\n"), "cov", 1);
 assert_checkerror("cov(""r"", [1;2], 1)", refMsg);
 
-assert_checkfalse(execstr("cov([1;2], ""r"", 1)"   ,"errcatch") == 0);
-refMsg = msprintf(_("%s: Wrong type for input argument #%d: a real matrix expected.\n"),"cov", 2);
+assert_checkfalse(execstr("cov([1;2], ""r"", 1)", "errcatch") =  = 0);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: a real matrix expected.\n"), "cov", 2);
 assert_checkerror("cov([1;2], ""r"", 1)", refMsg);
 
-assert_checkfalse(execstr("cov([1;2], [3;4], ""r"")"   ,"errcatch") == 0);
-refMsg = msprintf(_("%s: Wrong type for input argument #%d: an integer expected.\n"),"cov", 3);
+assert_checkfalse(execstr("cov([1;2], [3;4], ""r"")", "errcatch") =  = 0);
+refMsg = msprintf(_("%s: Wrong type for input argument #%d: an integer expected.\n"), "cov", 3);
 assert_checkerror("cov([1;2], [3;4], ""r"")", refMsg);
 
 x = [1;2];
 y = [3;4];
-computed = cov (x,y);
-expected = [0.5,0.5;0.5,0.5];
-assert_checkequal ( computed , expected );
+computed = cov(x, y);
+expected = [0.5, 0.5;0.5, 0.5];
+assert_checkalmostequal(computed, expected);
 //
-// The same, with nrmlztn=0
+// The same, with nrmlztn = 0
 x = [1;2];
 y = [3;4];
-computed = cov (x,y,0);
-expected = [0.5,0.5;0.5,0.5];
-assert_checkequal ( computed , expected );
+computed = cov(x, y, 0);
+expected = [0.5, 0.5;0.5, 0.5];
+assert_checkalmostequal(computed, expected);
 //
 x = [230;181;165;150;97;192;181;189;172;170];
 y = [125;99;97;115;120;100;80;90;95;125];
-expected = [1152.4556,-88.911111;-88.911111,244.26667];
-computed = cov (x,y);
-assert_checkalmostequal ( computed , expected , 1.e-7 );
+expected = [103721, -4001;-4001, 3664]./[90, 45;45, 15];
+computed = cov(x, y);
+assert_checkalmostequal(computed, expected);
 //
-// The same, with nrmlztn=0
+// The same, with nrmlztn = 0
 x = [230;181;165;150;97;192;181;189;172;170];
 y = [125;99;97;115;120;100;80;90;95;125];
-expected = [1152.4556,-88.911111;-88.911111,244.26667];
-computed = cov (x,y,0);
-assert_checkalmostequal ( computed , expected , 1.e-7 );
+expected = [103721, -4001;-4001, 3664]./[90, 45;45, 15];
+computed = cov(x, y, 0);
+assert_checkalmostequal(computed, expected);
 //
 x = [1;2;3;4;5];
-computed = cov (x);
+computed = cov(x);
 expected = 2.5;
-assert_checkequal ( computed , expected );
+assert_checkequal(computed, expected);
 //
-// The same, with nrmlztn=0
+// The same, with nrmlztn = 0
 x = [1;2;3;4;5];
-computed = cov (x,0);
+computed = cov(x, 0);
 expected = 2.5;
-assert_checkequal ( computed , expected );
+assert_checkequal(computed, expected);
 //
 A = [-1 1 2 ; -2 3 1 ; 4 0 3];
-Cexpected = [
-   10.3333   -4.1667    3.0000
-   -4.1667    2.3333   -1.5000
-    3.0000   -1.5000    1.0000
-];
-C = cov (A);
-assert_checkalmostequal ( Cexpected , C , [] , 1.e-4, "element");
+Cexpected =  [31, -25, 3;-25, 7, -3;3, -3, 1]./[3, 6, 1;6, 3, 2;1, 2, 1];
+C = cov(A);
+assert_checkalmostequal(Cexpected, C, [], [], "element");
 //
-// The same, with nrmlztn=0
+// The same, with nrmlztn = 0
 A = [-1 1 2 ; -2 3 1 ; 4 0 3];
-Cexpected = [
-   10.3333   -4.1667    3.0000
-   -4.1667    2.3333   -1.5000
-    3.0000   -1.5000    1.0000
-];
-C = cov(A,0);
-assert_checkalmostequal ( Cexpected , C , [] , 1.e-4, "element");
+C = cov(A, 0);
+assert_checkalmostequal(Cexpected, C, [], [], "element");
+//
+// Complex matrix
+A = complex([-1, 1, 2;-2, 3, 1;4, 0, 3], [2, -1, 1;1, -2, 3;3, 4, 0]);
+C = cov(A);
+// Diagonal has to be exactly real
+assert_checkequal(imag(diag(C)), [0;0;0])
+Cexpected = complex( [34, -7, 3;-7, 38, -17;3, -17, 10]./[3, 6, 2;6, 3, 3;2, 3, 3], ...
+                        [0, 71, -31;-71, 0, -2;31, 2, 0]./ [1, 6, 6;6, 1, 3;6, 3, 1])
+assert_checkalmostequal(C, Cexpected)
 //
 // Reference
 // 6.5.4.1. Mean Vector and Covariance Matrix
@@ -109,10 +110,10 @@ S = [
 0.0075 0.007 0.00135
 0.00175 0.00135 0.00043
 ];
-C = cov (A);
-assert_checkalmostequal ( S , C , 10*%eps , [] , "element");
+C = cov(A);
+assert_checkalmostequal(S, C, [], [], "element");
 //
-// The same, with nrmlztn=0
+// The same, with nrmlztn = 0
 A = [
 4.0 2.0 0.60
 4.2 2.1 0.59
@@ -125,31 +126,52 @@ S = [
 0.0075 0.007 0.00135
 0.00175 0.00135 0.00043
 ];
-C = cov (A,0);
-assert_checkalmostequal ( S , C , 10*%eps , [] , "element");
+C = cov(A, 0);
+assert_checkalmostequal(S, C, [], [], "element");
 //
 x = [1;2];
-computed = cov (x,1);
+computed = cov(x, 1);
 expected = 0.25;
-assert_checkequal ( computed , expected );
+assert_checkalmostequal(computed, expected);
 //
 x = [1;2];
-computed = cov (x,0);
+computed = cov(x, 0);
 expected = 0.5;
-assert_checkequal ( computed , expected );
+assert_checkalmostequal(computed, expected);
 //
 x = [1;2];
 y = [3;4];
-computed = cov (x,y,1);
-expected = [0.25,0.25;0.25,0.25];
-assert_checkequal ( computed , expected );
-//
-// Matlab compatibility
-x=[-1 1 2 ; -2 3 1 ; 4 0 3];
-computed = cov (x);
-expected = [
-   10.3333   -4.1667    3.0000
-   -4.1667    2.3333   -1.5000
-    3.0000   -1.5000    1.0000
-];
-assert_checkalmostequal ( computed , expected ,[],1.e-4);
+computed = cov(x, y, 1);
+expected = [0.25, 0.25;0.25, 0.25];
+assert_checkalmostequal(computed, expected);
+//
+// Case where x and y are matrices (cross-covariance)
+//
+x = [-1 1 2 ; -2 3 1 ; 4 0 3];
+y = [2, -1, 1;1, -2, 3;3, 4, 0]
+computed = cov(x, y);
+// check that cross-covariance is an exact submatrix of cov([x y])
+computed = cov(x, y);
+expected = cov([x, y])(1:3,4:$);
+assert_checkequal(computed, expected);
+// check actual value
+expected = [3,31,-25;-3,-25,7;1,3,-3]./[1,3,6;2,6,3;1,1,2];
+assert_checkalmostequal(computed, expected);
+computed = cov(x, y, 0);
+assert_checkalmostequal(computed, expected);
+computed = cov(x, y, 1);
+expected = [2,62,-25;-1,-25,14;2,2,-1]./[1,9,9;1,9,9;3,1,1];
+assert_checkalmostequal(computed, expected);
+//
+computed = cov(x, y);
+cov([x, y])(1:3,4:$)==cov(x,y)
+//
+assert_checkalmostequal(cov(1:5, 0), 2.5);
+assert_checkalmostequal(cov(1:5, 1), 2);
+assert_checkequal(cov(1:5, 0), cov((1:5)', 0));
+assert_checkequal(cov(1:5, 1), cov((1:5)', 1));
+//
+assert_checkequal(cov(2), 0);
+assert_checkequal(cov(2, 1), 0);
+assert_checkequal(cov(2, 2), zeros(2, 2));
+//