summaryrefslogtreecommitdiff
path: root/src/matrixOperations/division
diff options
context:
space:
mode:
Diffstat (limited to 'src/matrixOperations/division')
-rw-r--r--src/matrixOperations/division/Makefile.am18
-rw-r--r--src/matrixOperations/division/Makefile.in54
-rw-r--r--src/matrixOperations/division/crdiva.c33
-rw-r--r--src/matrixOperations/division/dldiva.c194
-rw-r--r--src/matrixOperations/division/drdiva.c268
-rw-r--r--src/matrixOperations/division/testMatrixDivision.c259
-rw-r--r--src/matrixOperations/division/testMatrixLDivision.c39
-rw-r--r--src/matrixOperations/division/testMatrixRDivision.c40
-rw-r--r--src/matrixOperations/division/zrdiva.c148
9 files changed, 443 insertions, 610 deletions
diff --git a/src/matrixOperations/division/Makefile.am b/src/matrixOperations/division/Makefile.am
index d17c95ec..37eb861f 100644
--- a/src/matrixOperations/division/Makefile.am
+++ b/src/matrixOperations/division/Makefile.am
@@ -11,7 +11,8 @@
##
libMatrixDivision_la_CFLAGS = -I $(top_builddir)/type \
- -I $(top_builddir)/matrixOperations/includes
+ -I $(top_builddir)/matrixOperations/includes \
+ -I $(top_builddir)/elementaryFunctions/includes
instdir = $(top_builddir)/lib
@@ -22,9 +23,11 @@ HEAD = ../includes/matrixDivision.h
libMatrixDivision_la_SOURCES = $(HEAD) \
srdiva.c \
drdiva.c \
- dldiva.c
+ dldiva.c \
+ zrdiva.c
-check_PROGRAMS = testMatrixLDivision
+check_PROGRAMS = testMatrixRDivision \
+ testMatrixLDivision
check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
$(top_builddir)/type/libFloatComplex.la \
@@ -36,8 +39,13 @@ check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
check_INCLUDES = -I $(top_builddir)/type \
-I $(top_builddir)/matrixOperations/includes
-testMatrixLDivision_SOURCES = testMatrixLDivision.c
+testMatrixLDivision_SOURCES = testMatrixLDivision.c
testMatrixLDivision_LDADD = $(check_LDADD)
testMatrixLDivision_CFLAGS = $(check_INCLUDES)
-TESTS = testMatrixLDivision
+testMatrixRDivision_SOURCES = testMatrixRDivision.c
+testMatrixRDivision_LDADD = $(check_LDADD)
+testMatrixRDivision_CFLAGS = $(check_INCLUDES)
+
+TESTS = testMatrixRDivision \
+ testMatrixLDivision
diff --git a/src/matrixOperations/division/Makefile.in b/src/matrixOperations/division/Makefile.in
index f3d456a9..eb093078 100644
--- a/src/matrixOperations/division/Makefile.in
+++ b/src/matrixOperations/division/Makefile.in
@@ -32,8 +32,9 @@ PRE_UNINSTALL = :
POST_UNINSTALL = :
build_triplet = @build@
host_triplet = @host@
-check_PROGRAMS = testMatrixLDivision$(EXEEXT)
-TESTS = testMatrixLDivision$(EXEEXT)
+check_PROGRAMS = testMatrixRDivision$(EXEEXT) \
+ testMatrixLDivision$(EXEEXT)
+TESTS = testMatrixRDivision$(EXEEXT) testMatrixLDivision$(EXEEXT)
subdir = matrixOperations/division
DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
@@ -56,7 +57,7 @@ libMatrixDivision_la_LIBADD =
am__objects_1 =
am_libMatrixDivision_la_OBJECTS = $(am__objects_1) \
libMatrixDivision_la-srdiva.lo libMatrixDivision_la-drdiva.lo \
- libMatrixDivision_la-dldiva.lo
+ libMatrixDivision_la-dldiva.lo libMatrixDivision_la-zrdiva.lo
libMatrixDivision_la_OBJECTS = $(am_libMatrixDivision_la_OBJECTS)
libMatrixDivision_la_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
$(LIBTOOLFLAGS) --mode=link $(CCLD) \
@@ -70,6 +71,14 @@ testMatrixLDivision_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
$(LIBTOOLFLAGS) --mode=link $(CCLD) \
$(testMatrixLDivision_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
$(LDFLAGS) -o $@
+am_testMatrixRDivision_OBJECTS = \
+ testMatrixRDivision-testMatrixRDivision.$(OBJEXT)
+testMatrixRDivision_OBJECTS = $(am_testMatrixRDivision_OBJECTS)
+testMatrixRDivision_DEPENDENCIES = $(check_LDADD)
+testMatrixRDivision_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
+ $(LIBTOOLFLAGS) --mode=link $(CCLD) \
+ $(testMatrixRDivision_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
+ $(LDFLAGS) -o $@
DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir)/includes
depcomp = $(SHELL) $(top_srcdir)/config/depcomp
am__depfiles_maybe = depfiles
@@ -83,9 +92,9 @@ LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \
--mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \
$(LDFLAGS) -o $@
SOURCES = $(libMatrixDivision_la_SOURCES) \
- $(testMatrixLDivision_SOURCES)
+ $(testMatrixLDivision_SOURCES) $(testMatrixRDivision_SOURCES)
DIST_SOURCES = $(libMatrixDivision_la_SOURCES) \
- $(testMatrixLDivision_SOURCES)
+ $(testMatrixLDivision_SOURCES) $(testMatrixRDivision_SOURCES)
ETAGS = etags
CTAGS = ctags
DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
@@ -200,7 +209,8 @@ target_alias = @target_alias@
top_builddir = @top_builddir@
top_srcdir = @top_srcdir@
libMatrixDivision_la_CFLAGS = -I $(top_builddir)/type \
- -I $(top_builddir)/matrixOperations/includes
+ -I $(top_builddir)/matrixOperations/includes \
+ -I $(top_builddir)/elementaryFunctions/includes
instdir = $(top_builddir)/lib
pkglib_LTLIBRARIES = libMatrixDivision.la
@@ -208,7 +218,8 @@ HEAD = ../includes/matrixDivision.h
libMatrixDivision_la_SOURCES = $(HEAD) \
srdiva.c \
drdiva.c \
- dldiva.c
+ dldiva.c \
+ zrdiva.c
check_LDADD = $(top_builddir)/type/libDoubleComplex.la \
$(top_builddir)/type/libFloatComplex.la \
@@ -223,6 +234,9 @@ check_INCLUDES = -I $(top_builddir)/type \
testMatrixLDivision_SOURCES = testMatrixLDivision.c
testMatrixLDivision_LDADD = $(check_LDADD)
testMatrixLDivision_CFLAGS = $(check_INCLUDES)
+testMatrixRDivision_SOURCES = testMatrixRDivision.c
+testMatrixRDivision_LDADD = $(check_LDADD)
+testMatrixRDivision_CFLAGS = $(check_INCLUDES)
all: all-am
.SUFFIXES:
@@ -295,6 +309,9 @@ clean-checkPROGRAMS:
testMatrixLDivision$(EXEEXT): $(testMatrixLDivision_OBJECTS) $(testMatrixLDivision_DEPENDENCIES)
@rm -f testMatrixLDivision$(EXEEXT)
$(testMatrixLDivision_LINK) $(testMatrixLDivision_OBJECTS) $(testMatrixLDivision_LDADD) $(LIBS)
+testMatrixRDivision$(EXEEXT): $(testMatrixRDivision_OBJECTS) $(testMatrixRDivision_DEPENDENCIES)
+ @rm -f testMatrixRDivision$(EXEEXT)
+ $(testMatrixRDivision_LINK) $(testMatrixRDivision_OBJECTS) $(testMatrixRDivision_LDADD) $(LIBS)
mostlyclean-compile:
-rm -f *.$(OBJEXT)
@@ -305,7 +322,9 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libMatrixDivision_la-dldiva.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libMatrixDivision_la-drdiva.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libMatrixDivision_la-srdiva.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libMatrixDivision_la-zrdiva.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testMatrixLDivision-testMatrixLDivision.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/testMatrixRDivision-testMatrixRDivision.Po@am__quote@
.c.o:
@am__fastdepCC_TRUE@ $(COMPILE) -MT $@ -MD -MP -MF $(DEPDIR)/$*.Tpo -c -o $@ $<
@@ -349,6 +368,13 @@ libMatrixDivision_la-dldiva.lo: dldiva.c
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libMatrixDivision_la_CFLAGS) $(CFLAGS) -c -o libMatrixDivision_la-dldiva.lo `test -f 'dldiva.c' || echo '$(srcdir)/'`dldiva.c
+libMatrixDivision_la-zrdiva.lo: zrdiva.c
+@am__fastdepCC_TRUE@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libMatrixDivision_la_CFLAGS) $(CFLAGS) -MT libMatrixDivision_la-zrdiva.lo -MD -MP -MF $(DEPDIR)/libMatrixDivision_la-zrdiva.Tpo -c -o libMatrixDivision_la-zrdiva.lo `test -f 'zrdiva.c' || echo '$(srcdir)/'`zrdiva.c
+@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/libMatrixDivision_la-zrdiva.Tpo $(DEPDIR)/libMatrixDivision_la-zrdiva.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='zrdiva.c' object='libMatrixDivision_la-zrdiva.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@ $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libMatrixDivision_la_CFLAGS) $(CFLAGS) -c -o libMatrixDivision_la-zrdiva.lo `test -f 'zrdiva.c' || echo '$(srcdir)/'`zrdiva.c
+
testMatrixLDivision-testMatrixLDivision.o: testMatrixLDivision.c
@am__fastdepCC_TRUE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixLDivision_CFLAGS) $(CFLAGS) -MT testMatrixLDivision-testMatrixLDivision.o -MD -MP -MF $(DEPDIR)/testMatrixLDivision-testMatrixLDivision.Tpo -c -o testMatrixLDivision-testMatrixLDivision.o `test -f 'testMatrixLDivision.c' || echo '$(srcdir)/'`testMatrixLDivision.c
@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/testMatrixLDivision-testMatrixLDivision.Tpo $(DEPDIR)/testMatrixLDivision-testMatrixLDivision.Po
@@ -363,6 +389,20 @@ testMatrixLDivision-testMatrixLDivision.obj: testMatrixLDivision.c
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixLDivision_CFLAGS) $(CFLAGS) -c -o testMatrixLDivision-testMatrixLDivision.obj `if test -f 'testMatrixLDivision.c'; then $(CYGPATH_W) 'testMatrixLDivision.c'; else $(CYGPATH_W) '$(srcdir)/testMatrixLDivision.c'; fi`
+testMatrixRDivision-testMatrixRDivision.o: testMatrixRDivision.c
+@am__fastdepCC_TRUE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixRDivision_CFLAGS) $(CFLAGS) -MT testMatrixRDivision-testMatrixRDivision.o -MD -MP -MF $(DEPDIR)/testMatrixRDivision-testMatrixRDivision.Tpo -c -o testMatrixRDivision-testMatrixRDivision.o `test -f 'testMatrixRDivision.c' || echo '$(srcdir)/'`testMatrixRDivision.c
+@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/testMatrixRDivision-testMatrixRDivision.Tpo $(DEPDIR)/testMatrixRDivision-testMatrixRDivision.Po
+@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='testMatrixRDivision.c' object='testMatrixRDivision-testMatrixRDivision.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixRDivision_CFLAGS) $(CFLAGS) -c -o testMatrixRDivision-testMatrixRDivision.o `test -f 'testMatrixRDivision.c' || echo '$(srcdir)/'`testMatrixRDivision.c
+
+testMatrixRDivision-testMatrixRDivision.obj: testMatrixRDivision.c
+@am__fastdepCC_TRUE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixRDivision_CFLAGS) $(CFLAGS) -MT testMatrixRDivision-testMatrixRDivision.obj -MD -MP -MF $(DEPDIR)/testMatrixRDivision-testMatrixRDivision.Tpo -c -o testMatrixRDivision-testMatrixRDivision.obj `if test -f 'testMatrixRDivision.c'; then $(CYGPATH_W) 'testMatrixRDivision.c'; else $(CYGPATH_W) '$(srcdir)/testMatrixRDivision.c'; fi`
+@am__fastdepCC_TRUE@ mv -f $(DEPDIR)/testMatrixRDivision-testMatrixRDivision.Tpo $(DEPDIR)/testMatrixRDivision-testMatrixRDivision.Po
+@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='testMatrixRDivision.c' object='testMatrixRDivision-testMatrixRDivision.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(testMatrixRDivision_CFLAGS) $(CFLAGS) -c -o testMatrixRDivision-testMatrixRDivision.obj `if test -f 'testMatrixRDivision.c'; then $(CYGPATH_W) 'testMatrixRDivision.c'; else $(CYGPATH_W) '$(srcdir)/testMatrixRDivision.c'; fi`
+
mostlyclean-libtool:
-rm -f *.lo
diff --git a/src/matrixOperations/division/crdiva.c b/src/matrixOperations/division/crdiva.c
index 02f8dc49..a3e528df 100644
--- a/src/matrixOperations/division/crdiva.c
+++ b/src/matrixOperations/division/crdiva.c
@@ -12,33 +12,10 @@
/**** WARNING NOT WORK AT ALL FOR THE MOMENT ***/
/*
-void crdiva ( floatComplex * in1, int it1,
- floatComplex * in2, int it2,
- floatComplex * out, int itOut, int size){
- int iIndex = 0; //Main loop index
- int iIndex1 = 0; //Loop index on left operand
- int iIndex2 = 0; //Loop index on right operand
- int iIndexOut = 0; //Lopp index on result matrix
+ void crdiva ( floatComplex * in1, int it1, floatComplex * in2, int it2, floatComplex * out) {
- for(iIndex = 0 ; iIndex < size ; iIndex++)
- {
- out[iIndexOut] = cdivides( in1[iIndex1], in2[iIndex2] ) ;
-
- iIndexOut += itOut;
- iIndex1 += it1;
- iIndex2 += it2;
- }
-
-
-}
-
-
-void cldiva ( floatComplex * in1, int it1,
- floatComplex * in2, int it2,
- floatComplex * out, int itOut, int size)) {
-f
-
- crdiva ( in2 , it2 , in1 , it1 , out , itout , size );
-}
-*/
+
+
+
+}*/
diff --git a/src/matrixOperations/division/dldiva.c b/src/matrixOperations/division/dldiva.c
index 09374d44..88b52f64 100644
--- a/src/matrixOperations/division/dldiva.c
+++ b/src/matrixOperations/division/dldiva.c
@@ -12,130 +12,104 @@
+
#include "matrixDivision.h"
#include "lapack.h"
#include <stdio.h>
#include <string.h>
-void dldiva ( double* in1, int lines1, int columns1 ,double* in2, int lines2, int columns2 ,double* out ){
-
- int iExit = 0 ;
- int info = 0 ;
- int iWork = 0 ;
- int tMinLinCol1 = 0 ;
- int tMaxLinCol1 = 0 ;
- int temp = 0 ;
- int rank = 0 ;
-
-
- int* pIpiv = NULL ;
- int* pIwork = NULL;
+void dldiva (double* in1, int lines1, int columns1 ,
+ double* in2, int lines2, int columns2 ,
+ double* out ){
+
+
+ char cNorm = 0;
+ int iExit = 0;
+
+ /*temporary variables*/
+ int iWork = 0;
+ int iInfo = 0;
+ int iMax = 0;
+ double dblRcond = 0;
+
+ double dblEps = 0;
+ double dblAnorm = 0;
+
+ double *pAf = NULL;
+ double *pXb = NULL;
+ double *pDwork = NULL;
+
+ int *pRank = NULL;
+ int *pIpiv = NULL;
int *pJpvt = NULL;
+ int *pIwork = NULL;
+
+ iWork = Max(4 * columns1, Max(Min(lines1, columns1) + 3 * lines1 + 1, 2 * Min(lines1, columns1) + columns2));
+
+
+
+
+ /* Array allocations*/
+ pAf = (double*)malloc(sizeof(double) * (unsigned int) lines1 * (unsigned int) columns1);
+ pXb = (double*)malloc(sizeof(double) * (unsigned int) Max(lines1,columns1) * (unsigned int) columns1);
+
+ pRank = (int*)malloc(sizeof(int));
+ pIpiv = (int*)malloc(sizeof(int) *(unsigned int) columns1);
+ pJpvt = (int*)malloc(sizeof(int) *(unsigned int) columns1);
+ pIwork = (int*)malloc(sizeof(int) *(unsigned int) columns1);
- char cNorm = 0;
-
- double rcond = 0;
- double anorm = 0;
- double epsilon = 0;
-
- double* cpyOfIn1 = NULL ;
- double* modifyIn2 = NULL ;
- double* pDwork = NULL ;
-
-
-/* adaptation of the original code in scilab
- iWorkMin = Max(4 * _iCols1, Max(Min(_iRows1, _iCols1) + 3 * _iRows1 + 1,
- 2 * Min(_iRows1, _iCols1) + _iRows2));
-*/
- if (lines1 > columns1)
- {
- tMinLinCol1 = columns1 ;
- tMaxLinCol1 = lines1 ;
- }
- else
- {
- tMinLinCol1 = lines1;
- tMaxLinCol1 = columns1;
- }
-
- if ( (tMinLinCol1 + 3 * lines1 + 1) > (2*tMinLinCol1 + lines2))
- temp = tMinLinCol1 + 3 * lines1 + 1;
- else
- temp = 2*tMinLinCol1 + lines2;
-
- if ( 4* columns1 > temp )
- iWork = 4* columns1 ;
- else
- iWork = temp ;
-
- pDwork = (double*)malloc(sizeof(double) * (unsigned int) iWork);
-
- cpyOfIn1 = (double*)malloc(sizeof(double) * (unsigned int) lines1
- * (unsigned int) columns1);
- modifyIn2 = (double*)malloc(sizeof(double) * (unsigned int) tMaxLinCol1
- * (unsigned int) columns1);
-
- pIpiv = (int*)malloc(sizeof(int) * (unsigned int) columns1);
- pIwork = (int*)malloc(sizeof(int) * (unsigned int) columns1);
- pJpvt = (int*)malloc(sizeof(int) * (unsigned int) columns1);
-
- epsilon = getRelativeMachinePrecision() ;
-
+
+
cNorm = '1';
- anorm = dlange_ (&cNorm , &lines1, &columns1, in1, &lines1, pDwork);
+ pDwork = (double*)malloc(sizeof(double) *(unsigned int)iWork);
+ dblEps = getRelativeMachinePrecision() ;
- if ( lines1 == columns1 )
+ dblAnorm = dlange_(&cNorm, &lines1, &columns1, in1, &lines1, pDwork);
+
+ if(lines1 == columns1)
+ {
+ cNorm = 'F';
+ dlacpy_(&cNorm, &columns1, &columns1, in1, &columns1, pAf, &columns1);
+ dgetrf_(&columns1, &columns1, pAf, &columns1, pIpiv, &iInfo);
+ if(iInfo == 0)
{
- cNorm = 'F' ;
- dlacpy_ ( &cNorm , &columns1 , &columns1, in1, &columns1, cpyOfIn1,
- &columns1 );
- dgetrf_ (&columns1, &columns1, cpyOfIn1, &columns1, pIpiv, &info);
-
- if ( info == 0 )
- {
- cNorm = '1' ;
- dgecon_ ( &cNorm, &columns1, cpyOfIn1, &columns1, &anorm,
- &rcond, pDwork, pIwork, &info);
-
- if ( rcond > sqrt (epsilon ))
- {
- cNorm = 'N' ;
- dgetrs_ ( &cNorm, &columns1, &columns2, cpyOfIn1, &columns1,
- pIpiv, in2, &columns1, &info);
- cNorm = 'F' ;
- dlacpy_ ( &cNorm, &columns1, &columns2, in2, &columns1,
- out, &columns1 );
- iExit = 1 ;
-
- }
- }
+ cNorm = '1';
+ C2F(dgecon)(&cNorm, &columns1, pAf, &columns1, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo);
+ if(dblRcond > sqrt(dblEps))
+ {
+ cNorm = 'N';
+ C2F(dgetrs)(&cNorm, &columns1, &columns2, pAf, &columns1, pIpiv, in2, &columns1, &iInfo);
+ cNorm = 'F';
+ C2F(dlacpy)(&cNorm, &columns1, &columns2, in2, &columns1, out, &columns1);
+ iExit = 1;
+ }
}
+ }
+ if(iExit == 0)
+ {
+ dblRcond = sqrt(dblEps);
+ cNorm = 'F';
+ iMax = Max(lines1, columns1);
+ C2F(dlacpy)(&cNorm, &lines1, &columns2, in2, &lines2, pXb, &iMax);
+ memset(pJpvt, 0x00,(unsigned int) sizeof(int) * (unsigned int) columns1);
+ C2F(dgelsy)( &lines1, &columns1, &columns2, in1, &lines1, pXb, &iMax,
+ pJpvt, &dblRcond, &pRank[0], pDwork, &iWork, &iInfo);
- if ( iExit == 0 )
+ if(iInfo == 0)
{
- rcond = sqrt ( epsilon ) ;
- cNorm = 'F' ;
- dlacpy_ ( &cNorm , &lines1, &columns2, in2, &lines1, modifyIn2, &tMaxLinCol1);
- memset ( pJpvt , 0x00, sizeof(int) * (unsigned int) columns1) ;
- dgelsy_ ( &lines1, &columns1, &columns2, in1, &lines1 ,
- modifyIn2 , &tMaxLinCol1, pJpvt, &rcond, &rank,
- pDwork, &iWork, &info );
-
- if (info ==0 )
- {
- cNorm = 'F' ;
- dlacpy_ ( &cNorm, &columns1, &columns2, modifyIn2, &tMaxLinCol1,
- out , &columns1 );
-
- }
+
+ cNorm = 'F';
+ C2F(dlacpy)(&cNorm, &columns1, &columns2, pXb, &iMax, out, &columns1);
}
-
- free (pIpiv) ;
- free (pIwork);
- free (pJpvt);
- free (cpyOfIn1 ) ;
- free (modifyIn2 ) ;
- free ( pDwork ) ;
+ }
+
+ free(pAf);
+ free(pXb);
+ free(pRank);
+ free(pIpiv);
+ free(pJpvt);
+ free(pIwork);
+ free(pDwork);
}
diff --git a/src/matrixOperations/division/drdiva.c b/src/matrixOperations/division/drdiva.c
index ab2318aa..bf94533d 100644
--- a/src/matrixOperations/division/drdiva.c
+++ b/src/matrixOperations/division/drdiva.c
@@ -10,164 +10,136 @@
*
*/
-
-
#include "matrixDivision.h"
#include "lapack.h"
+#include <string.h>
+#include <stdio.h>
+
+
+
+void drdiva ( double * in1, int lines1, int columns1,
+ double * in2, int lines2, int columns2,
+ double * out){
+
+ char cNorm = 0;
+ int iExit = 0;
+
+ /*temporary variables*/
+ int iWork = 0;
+ int iInfo = 0;
+ int iMax = 0;
+ double dblRcond = 0;
+
+ double dblEps = 0;
+ double dblAnorm = 0;
+
+ double *pAf = NULL;
+ double *pAt = NULL;
+ double *pBt = NULL;
+ double *pDwork = NULL;
+
+ int *pRank = NULL;
+ int *pIpiv = NULL;
+ int *pJpvt = NULL;
+ int *pIwork = NULL;
+
+ iWork = Max(4 * columns2, Max(Min(lines2, columns2) + 3 * lines2 + 1, 2 * Min(lines2, columns2) + lines1));
+
+
+ /* Array allocations*/
+ pAf = (double*)malloc(sizeof(double) * (unsigned int)columns2 * (unsigned int)lines2);
+ pAt = (double*)malloc(sizeof(double) * (unsigned int)columns2 *(unsigned int) lines2);
+ pBt = (double*)malloc(sizeof(double) * (unsigned int)Max(lines2,columns2) * (unsigned int)lines1);
+
+ pRank = (int*)malloc(sizeof(int));
+ pIpiv = (int*)malloc(sizeof(int) * (unsigned int)columns2);
+ pJpvt = (int*)malloc(sizeof(int) * (unsigned int)lines2);
+ pIwork = (int*)malloc(sizeof(int) * (unsigned int)columns2);
+
+
+ cNorm = '1';
+ pDwork = (double*)malloc(sizeof(double) * (unsigned int)iWork);
+ dblEps = getRelativeMachinePrecision() ;
+ dblAnorm = dlange_(&cNorm, &lines2, &columns1, in2, &lines2, pDwork);
-void drdiva ( double* in1, int lines1, int columns1 ,
- double* in2, int lines2, int columns2 ,
- double* out ){
-
-/* creation of all temp variable , maybe some of them are not accurate and could be deleted
-**/
-
- int iexit = 0 ;
- int iwork = 0 ;
- int info = 0 ;
- int imax = 0 ;
- int rank = 0 ;
-
- int tMinLinCol1 = 0 ;
- int temp = 0;
-
- int i = 0 , j = 0 ,ij = 0 ,ji = 0;
-
-
- int* pIpiv = (int*) malloc( sizeof(int) * (unsigned int) columns1);
- int* pIwork = (int*) malloc( sizeof(int) * (unsigned int) columns1);
- int* pJpvt = (int*) malloc( sizeof(int) * (unsigned int) lines1);
-
- double* work = NULL;
- double anorm = 0;
- double rcond = 0 ;
- double epsilon = getRelativeMachinePrecision() ;
-
-
-
-
- double* transpOfIn1 = (double*) malloc(sizeof(double) * (unsigned int) lines1
- * (unsigned int) columns1);
-
- double* transpOfIn2 = (double*) malloc(sizeof(double) * (unsigned int) lines2
- * (unsigned int) columns2);
-
- double* copyOfTransIn1 = (double*) malloc( sizeof(double) * (unsigned int) lines1
- * (unsigned int) columns1);
-
-
-
-/* adaptation of the original code in scilab
- iWorkMin = Max(4 * _iCols1, Max(Min(_iRows1, _iCols1) + 3 * _iRows1 + 1,
- 2 * Min(_iRows1, _iCols1) + _iRows2));
- */
- if (lines1 > columns1)
- tMinLinCol1 = columns1 ;
- else
- tMinLinCol1 = lines1;
-
- if ( (tMinLinCol1 + 3 * lines1 + 1) > (2*tMinLinCol1 + lines2))
- temp = tMinLinCol1 + 3 * lines1 + 1;
- else
- temp = 2*tMinLinCol1 + lines2;
-
- if ( 4* columns1 > temp )
- iwork = 4* columns1 ;
- else
- iwork = temp ;
-
-
-
-
- work = (double*) malloc (sizeof(double) *(unsigned int) iwork );
-
- anorm = getOneNorm( &lines1, &columns1 , in1 , work);
-
-
-/* end of allocation area*/
-
-
- dtransposea ( in1, lines1, columns1, transpOfIn1);
- dtransposea ( in2, lines2, columns2, transpOfIn2) ;
-
-
-
-/* case of a square matrix */
- if ( lines1 == columns1 )
- {
- drowcata ( transpOfIn1, lines1, columns1, NULL, 0 , 0 , copyOfTransIn1 ) ;
-
-
- /*/ put here algo of LU fact of in1
- dgetrf ( &columns1 , &columns1 , in1 , &columns1 , pIpiv , &info )
- //return value in pIpiv*/
- dgetrf_ ( &columns1 , &columns1 ,copyOfTransIn1 , &columns1 , pIpiv , &info );
-
- if ( info == 0 )
+ /*tranpose A and B*/
+
+ dtransposea(in2, lines2, columns2, pAt);
+ dtransposea(in1, lines1, columns2, pBt);
+
+ if(lines2 == columns2)
+ {
+ cNorm = 'F';
+ dlacpy_(&cNorm, &columns2, &columns2, pAt, &columns2, pAf, &columns2);
+ dgetrf_(&columns2, &columns2, pAf, &columns2, pIpiv, &iInfo);
+ if(iInfo == 0)
{
- /*/to get 1-norm of in1 put here algo dgecon which return rcond*/
- dgecon_("1" , &columns1, copyOfTransIn1 , &columns1, &anorm,
- &rcond, work, pIwork, &info);
-
- if ( rcond > sqrt(epsilon ))
- {/* put here algo to resolv linear equation in1 * X = in2 ,
- the return value go in in2
- // put here algo to copy in2 in out */
- resolveSystemLinear (&columns1, &lines2, copyOfTransIn1, pIpiv,
- transpOfIn2, &info) ;
- dtransposea ( transpOfIn2 , columns1 , lines2 , out) ;
- iexit = 1;
- }
-
+ cNorm = '1';
+ dgecon_(&cNorm, &columns2, pAf, &columns2, &dblAnorm, &dblRcond, pDwork, pIwork, &iInfo);
+ if(dblRcond > sqrt(dblEps))
+ {
+ cNorm = 'N';
+ dgetrs_(&cNorm, &columns2, &lines1, pAf, &columns2, pIpiv, pBt, &columns2, &iInfo);
+ dtransposea(pBt, columns2, lines1, out);
+ iExit = 1;
+ }
}
-
}
-
-/* non-square matrix case */
- if ( iexit== 0)
+ if(iExit == 0)
{
- rcond = sqrt(epsilon );
- if (lines1 < columns1)
- imax = columns1 ;
- else
- imax = lines1;
-
-
-
- dgelsy_ (&columns1, &lines1, &lines2, transpOfIn1 , &columns1,
- transpOfIn2, &imax, pJpvt, &rcond, &rank, work,
- &iwork, &info);
-
-
- if (info == 0)
- for(j = 0 ; j < lines1 ; j++)
+ dblRcond = sqrt(dblEps);
+ cNorm = 'F';
+ iMax = Max(lines2, columns2);
+ memset(pJpvt, 0x00, (unsigned int)sizeof(int) * (unsigned int)lines2);
+ dgelsy_(&columns2, &lines2, &lines1, pAt, &columns2, pBt, &iMax,
+ pJpvt, &dblRcond, &pRank[0], pDwork, &iWork, &iInfo);
+
+ if(iInfo == 0)
+ {
+
+
+ /* TransposeRealMatrix(pBt, lines1, lines2, out, Max(lines1,columns1), lines2);*/
+
+ /*Mega caca de la mort qui tue des ours a mains nues
+ mais je ne sais pas comment le rendre "beau" :(*/
+ {
+ int i,j,ij,ji;
+ for(j = 0 ; j < lines2 ; j++)
{
- for(i = 0 ; i < lines2 ; i++)
- {
- ij = i + j * lines2;
- if ( lines1 > columns1 )
- ji = j + i *lines1 ;
- else
- ji = j + i *columns1;
- out[ij] = transpOfIn2[ji];
- }
+ for(i = 0 ; i < lines1 ; i++)
+ {
+ ij = i + j * lines1;
+ ji = j + i * Max(lines2, columns2);
+ out[ij] = pBt[ji];
+ }
}
-
-
+ }
+ }
}
-
-
-
- /* then we free all the allocated memory */
-
- free ( pIpiv ) ;
- free ( pIwork) ;
- free ( work ) ;
-
- free (transpOfIn1) ;
- free (transpOfIn2) ;
- free (copyOfTransIn1);
-
+
+ free(pAf);
+ free(pAt);
+ free(pBt);
+ free(pRank);
+ free(pIpiv);
+ free(pJpvt);
+ free(pIwork);
+ free(pDwork);
+
+}
+
+int Max(int _dblVar1, int _dblVar2)
+{
+ if(_dblVar1 > _dblVar2)
+ return _dblVar1;
+ return _dblVar2;
+}
+
+int Min(int _dblVar1, int _dblVar2)
+{
+ if(_dblVar1 < _dblVar2)
+ return _dblVar1;
+ return _dblVar2;
}
diff --git a/src/matrixOperations/division/testMatrixDivision.c b/src/matrixOperations/division/testMatrixDivision.c
deleted file mode 100644
index ba3e1839..00000000
--- a/src/matrixOperations/division/testMatrixDivision.c
+++ /dev/null
@@ -1,259 +0,0 @@
-/*
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008-2008 - INRIA - Allan SIMON
- *
- * This file must be used under the terms of the CeCILL.
- * This source file is licensed as described in the file COPYING, which
- * you should have received as part of this distribution. The terms
- * are also available at
- * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
- *
- */
-
-#include "matrixDivision.h"
-#include <assert.h>
-#include <stdio.h>
-
-#define LINES1 2
-#define LINES2 2
-#define COLUMNS 2
-/*
-static void sdivaTest ( void )
-{
- int i = 0 ;
-
- float in1[] =
- {0.51323400903493166f,0.6807207581587136f,0.54669387824833393f,0.24904179340228438f,
- 0.95017496403306723f,0.91187966475263238f,0.78931419923901558f,0.30985609395429492f,
- 0.71984737459570169f,0.9819172923453152f,0.76060794852674007f,0.504620002117008f,
- 0.06844846438616514f,0.1693508871831f,0.76318027824163437f,0.30514528928324580f,
- 0.92679917532950640f,0.04813073994591832f,0.93056132830679417f,0.31760499393567443f,
-0.20109100220724940f,0.75775502342730761f,0.15586951794102788f,0.59756303764879704f,
- 0.93055095290765166f,0.94091763999313116f,0.42790159443393350f,0.01439402624964714f,
- 0.85971397114917636f,0.11901073250919580f,0.59864782588556409f,0.15444914810359478f,
- 0.44080717256292701f,0.09636751096695662f,0.47461007768288255f,0.42803008854389191f,
- 0.92924218205735087f,0.47864412236958742f,0.44434435339644551f,0.52509398944675922f,
-0.78608208894729614f,0.46497652260586619f,0.06789979804307222f,0.24485790403559804f,
- 0.71605867333710194f,0.99458031123504043f,0.84305586572736502f,0.45727639505639672f,
- 0.29075053706765175f,0.55482550663873553f,0.28596154693514109f,0.07587631093338132f,
- 0.66903266869485378f,0.32733985921368f,0.90166416298598051f,0.83553476119413972f,
- 0.80835641175508499f,0.73590047238394618f,0.83198319096118212f,0.93355408729985356f,
-0.78650354826822877f,0.74095427244901657f,0.95994638977572322f,0.51785656530410051f,
- 0.78444739105179906f,0.13383972086012363f,0.43243235861882567f,0.74689115490764380f,
- 0.96420694747939706f,0.22162469848990440f,0.41627834690734744f,0.81930211279541254f,
- 0.21618453459814191f,0.08386900834739208f,0.48852836480364203f,0.20842899661511183f,
- 0.22902313107624650f,0.78962677717208862f,0.25109061924740672f,0.57858273852616549f,
-0.69515300076454878f,0.48902340466156602f,0.47249071300029755f,0.59868981270119548f,
- 0.75798543263226748f,0.80475882859900594f,0.61186199076473713f,0.6943939602933824f,
- 0.10014689620584249f,0.01250550756230950f,0.47080435231328011f,0.58039451343938708f,
- 0.95630660001188517f,0.32801365898922086f,0.27026554010808468f,0.52016736706718802f,
- 0.16067446302622557f,0.04408275568857789f,0.88035558909177780f,0.92851745663210750f,
-0.42481321236118674f,0.73924016486853361f,0.16896168375387788f,0.39154956489801407f,
- 0.97163037536665797f,0.88981838244944811f,0.90427244128659368f,0.31537816859781742f,
- 0.06473635649308562f,0.58296835329383612f,0.14059370616450906f,0.63762421533465385f,
- 0.67373040271922946f,0.44792105350643396f,0.60066422121599317f,0.06632651202380657f,
- 0.66453591873869300f,0.19733488839119673f,0.45683057839050889f,0.08710412681102753f,
-0.34466254524886608f,0.09304937114939094f,0.06055234652012587f,0.10790407890453935f,
- 0.15703585743904114f,0.72192603675648570f,0.34105927217751741f,0.62482782872393727f,
- 0.68893781490623951f,0.03671516245231032f,0.98292266484349966f,0.73650254914537072f,
- 0.13196587935090065f,0.18807678623124957f,0.43376339320093393f,0.47929613338783383f,
- 0.15195304714143276f,0.18539744755253196f,0.92726647388190031f,0.80549291754141450f,
-0.64977857517078519f,0.70831089280545712f,0.90016864379867911f,0.51403949689120054f,
- 0.49954565847292542f,0.55090149492025375f,0.92046913085505366f,0.74058383423835039f,
- 0.82662396552041173f,0.3313873652368784f,0.68757036840543151f,0.06999884452670813f,
- 0.41287241736426950f,0.04925781115889549f,0.85545881045982242f,0.31712341401726007f,
- 0.03992868261411786f,0.92384314350783825f,0.74299975624307990f,0.00424729567021132f,
-0.01225362811237574f,0.3223448325879872f,0.93296645395457745f,0.08063758304342628f,
- 0.74957344215363264f,0.82469086581841111f,0.35314525663852692f,0.44088636664673686f,
- 0.79343967605382204f,0.06321920128539205f,0.87087013013660908f,0.05352633958682418f,
- 0.47160778101533651f,0.38195306668058038f,0.04547535255551338f,0.02313599688932300f,
- 0.63996278587728739f,0.36444053100422025f,0.37049167416989803f,0.76907502254471183f,
-0.89965870184823871f,0.89933154825121164f,0.64569224463775754f,0.34696785174310207f,
- 0.39040711661800742f,0.08694788347929716f,0.22625351930037141f,0.21681279689073563f,
- 0.33453882811591029f,0.15584628004580736f,0.95653126062825322f,0.83982629515230656f,
- 0.50479181623086333f,0.54758223798125982f,0.83060362795367837f,0.21214072033762932f,
- 0.02860224200412631f,0.95791505370289087f,0.91556971566751599f,0.94719038717448711f,
-0.30791273340582848f,0.81771020544692874f,0.74679336044937372f,0.87552759842947125f,
- 0.49545058421790600f,0.48191254725679755f,0.39087839704006910f,0.88425681227818131f,
- 0.54299664497375488f,0.20552197424694896f,0.31019500363618135f,0.57183724315837026f,
- 0.54999292083084583f,0.12055991357192397f,0.26517685409635305f,0.29529260704293847f,
- 0.95132300630211830f,0.57458581728860736f,0.71528563741594553f,0.91191364871338010f,
-0.18359116325154901f,0.45608301833271980f,0.17441136343404651f,0.18253823462873697f,
- 0.77341705607250333f,0.01534702442586422f,0.27907355269417167f,0.94600243400782347f,
- 0.39083331311121583f,0.35583620518445969f,0.58566563902422786f,0.44004907924681902f,
- 0.79277362348511815f,0.96740394271910191f,0.06400812184438109f,0.07406814303249121f,
- 0.03709788480773568f,0.85063817724585533f,0.10404936922714114f,0.12005183193832636f,
-0.52963322307914495f,0.42104291776195168f,0.92252827808260918f,0.17444357229396701f,
- 0.34813721571117640f,0.98178615467622876f,0.51038642041385174f,0.35229418566450477f,
- 0.33187932055443525f,0.12522496515884995f,0.17751775681972504f,0.30953403143212199f,
- 0.30514361429959536f,0.69633625121787190f,0.58739017136394978f,0.95729830628260970f,
- 0.78135449346154928f,0.41670671710744500f,0.55824907496571541f,0.17518991930410266f,
-0.11329598492011428f,0.87877958174794912f,0.81112976977601647f,0.1327551044523716f,
- 0.88772260351106524f,0.77383322361856699f,0.39891980635002255f,0.42609489522874355f,
- 0.29346287390217185f,0.29799025785177946f,0.62878308678045869f,0.90327445417642593f,
- 0.88942573545500636f,0.11638559121638536f,0.92788035096600652f,0.30957929231226444f,
- 0.2565767071209848f,0.42655616905540228f,0.26941573480144143f,0.65786541625857353f,
-0.98608913458883762f,0.04273471748456359f,0.26294819917529821f,0.66961710015311837f,
- 0.19781696423888206f,0.62439860356971622f,0.25354105327278376f,0.55396229820325971f,
- 0.92479544691741467f,0.09095242014154792f,0.63873832207173109f,0.97996837133541703f,
- 0.67638632655143738f,0.89037371007725596f,0.41085386741906404f,0.07612052233889699f,
- 0.31289586611092091f,0.62697393959388137f,0.15661530848592520f,0.28013094374909997f};
-
-float in2[] =
- {0.56034345272928476f,0.89566554129123688f,0.53930272068828344f,0.3547350224107504f,
- 0.80800013709813356f,0.62323769554495811f,0.58200186025351286f,0.56180121190845966f,
- 0.23549679014831781f,0.28873602300882339f,0.43252215441316366f,0.7673156540840864f,
- 0.06873596925288439f,0.18797885254025459f,0.86748637538403273f,0.73920361138880253f,
- 0.97392784897238016f,0.86080306768417358f,0.39093428757041693f,0.57339327596127987f,
-0.16198171628639102f,0.66783405328169465f,0.1454864419065416f,0.32747871475294232f,
- 0.51621831534430385f,0.79978153714910150f,0.09275748720392585f,0.14153907122090459f,
- 0.06119967205449939f,0.32070356840267777f,0.73968251561746001f,0.37837028549984097f,
- 0.56752133695408702f,0.11224916437640786f,0.28770424565300345f,0.00037088664248586f,
- 0.79149663401767612f,0.58377730334177613f,0.59390504425391555f,0.02694623963907361f,
-0.76592414453625679f,0.02564378362149000f,0.74512455798685551f,0.02366107050329447f,
- 0.06156063079833984f,0.04258572962135077f,0.63941287063062191f,0.38405111897736788f,
- 0.14321060106158257f,0.42103306483477354f,0.05339348502457142f,0.93415357265621424f,
- 0.59929492324590683f,0.66478141304105520f,0.14329732768237591f,0.42867958266288042f,
- 0.10990926995873451f,0.63394964020699263f,0.05351450480520725f,0.54702291730791330f,
-0.69746216991916299f,0.03159578284248710f,0.25769635709002614f,0.59392183972522616f,
- 0.01609914982691407f,0.86753786867484450f,0.32879876391962171f,0.22861831961199641f,
- 0.33929981896653771f,0.75886590173467994f,0.61312689306214452f,0.48855357570573688f,
- 0.10645245248451829f,0.14542592084035277f,0.18745915638282895f,0.87820987729355693f,
- 0.92918653646484017f,0.92348486324772239f,0.39261205168440938f,0.68466226710006595f,
-0.08946218248456717f,0.42732305638492107f,0.77090662438422441f,0.924068246036768f,
- 0.50982708018273115f,0.34508761204779148f,0.10317245963960886f,0.63696919381618500f,
- 0.88041578140109777f,0.65034613572061062f,0.80938913393765688f,0.9986613355576992f,
- 0.38141551148146391f,0.43060396797955036f,0.61948752496391535f,0.80996788293123245f,
- 0.05324298795312643f,0.29668187908828259f,0.00448737759143114f,0.82276185229420662f,
-0.19805425917729735f,0.10866974340751767f,0.99788628844544291f,0.6723356381990016f,
- 0.46966064115986228f,0.79823006363585591f,0.67132972134277225f,0.35570297623053193f,
- 0.52579802041873336f,0.25843874411657453f,0.19287035940214992f,0.59324032673612237f,
- 0.25955950608476996f,0.61404782952740788f,0.21750316722318530f,0.81994143361225724f,
- 0.19624035572633147f,0.8358787004835904f,0.42109713284298778f,0.31491625169292092f,
-0.25988535769283772f,0.76794129703193903f,0.84549946337938309f,0.48446214850991964f,
- 0.52822500281035900f,0.37241784948855639f,0.85067357495427132f,0.16846220474690199f,
- 0.71482414938509464f,0.78300847951322794f,0.54156896471977234f,0.89239248540252447f,
- 0.00643130205571651f,0.21658254135400057f,0.44520513340830803f,0.32618630956858397f,
- 0.52326664514839649f,0.32250450644642115f,0.26230763643980026f,0.23438148852437735f,
-0.50999558391049504f,0.37397424085065722f,0.64668390387669206f,0.64172910666093230f,
- 0.03306737588718534f,0.18668571440503001f,0.83025926211848855f,0.99121205648407340f,
- 0.38945918949320912f,0.27384403301402926f,0.41589357936754823f,0.99599931901320815f,
- 0.81049045221880078f,0.96830060658976436f,0.40803860733285546f,0.52525822212919593f,
- 0.33058117749169469f,0.23706211848184466f,0.17479355866089463f,0.62723324215039611f,
-0.16741782892495394f,0.15187738463282585f,0.72003478836268187f,0.59474316425621510f,
- 0.71126131806522608f,0.50548844784498215f,0.43874060269445181f,0.84671537391841412f,
- 0.46396317798644304f,0.84823036566376686f,0.28646126668900251f,0.48104315437376499f,
- 0.26604998949915171f,0.66378767788410187f,0.63686545100063086f,0.79531485401093960f,
- 0.95694970060139894f,0.70720722898840904f,0.83713256847113371f,0.62996550090610981f,
-0.34009417472407222f,0.72199993440881371f,0.67456434061750770f,0.10675506712868810f,
- 0.79017778439447284f,0.34282173449173570f,0.58765271818265319f,0.07389529095962644f,
- 0.16408033994957805f,0.12704358855262399f,0.51104495069012046f,0.07703803153708577f,
- 0.12485344661399722f,0.51966900611296296f,0.49762418633326888f,0.77279568510130048f,
- 0.95798523304983974f,0.22388020763173699f,0.83013197174295783f,0.04275623383000493f};
-
-
-float result[] =
-{-0.13425098694507090f,0.24735018618193444f,-0.13997367470679392f,0.11987585715392905f,
- -0.07023550140255914f,-0.32818483956327921f,0.82293845909628660f,-0.04825208923032670f,
- 0.41107043821212774f,0.28834120500378574f,
-0.26097631689267364f,0.41112922481604613f,-0.48704150701173060f,0.16756167892414930f,
- 0.18525816194433214f,0.24194769618197037f,0.04516579999951376f,0.08210988946198612f,
- -0.17724539433971653f,0.21570365259908197f,
-0.46476238563747729f,0.12929328091152756f,0.03791742144738922f,0.33232469296313982f,
- -0.14531956716307623f,-0.04355937726216118f,-0.32871318476723438f,0.16208735115250936f,
- 0.27759742819012267f,0.20883772519304408f,
-0.60589399320385928f,0.06028369176400562f,0.33610957414517950f,-0.32741653379843827f,
- 0.32006732568590329f,0.12357809346654614f,0.26645037101516444f,-0.15778041923768410f,
- -0.02281069075548223f,-0.21831876232917027f,
-0.11316162639341981f,0.44869069077875356f,0.05087241573803265f,-0.07310286041500898f,
- 0.24172631454502305f,0.20582249789901347f,-0.49030030296647070f,0.49725834464740676f,
- -0.15098926456955386f,0.29845287904710704f,
-0.26126025652717821f,0.61602266860724408f,-0.04395376298763495f,-0.31468851985942192f,
- -0.42610430769336133f,0.07619666627024707f,0.42141080474076514f,0.34732212135875368f,
- 0.21966932936729675f,-0.16465486570712187f,
-0.01801836442044520f,0.13867951307716378f,-0.22823955222461578f,0.35528257041935002f,
- 0.3416886047495104f,-0.05324845913867078f,-0.18334602929956423f,0.34486570917896908f,
- -0.09850219178062529f,0.15717539967151753f,
-0.26027320955453370f,0.24744146889140337f,-0.164185466873868f,-0.18897859269730288f,
- -0.10286928200189474f,0.54800666111794050f,0.32097861553203794f,0.16617893447929227f,
- -0.03433269255504845f,-0.01421332288783149f,
--0.24887241840873561f,0.52303584234440914f,0.02976416459895397f,-0.12264750867737524f,
- 0.30087087321455463f,0.28120264259279731f,-0.02981261204955734f,-0.05832040687613257f,
- 0.09480786369785668f,0.26640553810556677f,
-0.28112837735944118f,0.53844297257991636f,0.87478164836631844f,-0.1409682347398108f,
- 0.66145942836153393f,-0.26119505204825816f,-0.34297154024663806f,-0.24798246038042729f,
- -0.30432340939641339f,0.37887301625026215f,
-0.38327878821721156f,0.28922881585078497f,-0.34742629642806333f,-0.03677194909448467f,
- -0.04829379940842784f,0.06878313210496664f,0.04495930782786764f,0.47892284040726463f,
- 0.32833980213443120f,-0.23534654463679133f,
--0.02763634740677747f,0.72078725530634646f,0.48763002329592203f,-0.57876872035917137f,
- 0.23063708055466298f,0.04745273757978258f,-0.28751023415434818f,0.138229588096174f,
- 0.36392895348902060f,-0.15741181362164816f,
-0.13480297497174423f,-0.06082319210597478f,-0.03346206227374363f,0.157369564795072f,
- -0.09939329142432089f,0.50801256096451919f,-0.21552055455365590f,0.18726599728475304f,
- -0.09981549922929681f,0.56837678961938298f,
-0.51101691288362061f,0.72117205689181529f,0.47311124338037619f,-0.59816447223287172f,
- 0.56148184747810326f,0.06803537921378118f,-0.12236345370998958f,-0.2366952577209536f,
- -0.10560524465448745f,-0.03700000431544657f,
--0.07318875604760167f,0.49681103126620946f,0.16658812117729807f,0.19247716660649325f,
- 0.17479922937365650f,0.20305642078183334f,-0.11893261285117168f,0.52167729786941019f,
- -0.26590601210923126f,-0.17052720242327102f};
-
-float out[(COLUMNS)*(LINES2)] ;
-
- printf("\n>rentre dans fonction \n");
- srdiva ( in1 , LINES1 , COLUMNS , in2 , LINES2 , COLUMNS , out ) ;
- printf("\n\n\t>>>>>debut assert\n");
- for ( i = 0 ; i < LINES2 *COLUMNS ; i++ )
- {
- printf ( "\n %d out : %e result : %e\n" , i , out[i] , result[i] ) ;
- assert ( fabs ( out[i] - result[i] ) / fabs( out[i]) < 1e-6 ) ;
- }
-
-}
-
-
-
-
-*/
-
-static void ddivaTest ( void )
-{
- int i = 0 ;
-
- double in1[] = { 1 , 2 , 3 , 4 } ;
-
-
-double in2[] = { 1 , 2 ,3 ,4 } ;
-
-
-
-double result[] = { 1 ,0 , 0 , 1 };
-
-
-double out [(COLUMNS)*(LINES2)] ;
-
- drdiva ( in1 , LINES1 , COLUMNS , in2 , LINES2 , COLUMNS , out ) ;
- printf("\n\n\t>>>>>debut assert\n");
- for ( i = 0 ; i < LINES2 *COLUMNS ; i++ )
- {
- printf ( "\n %d out : %e result : %e\n" , i , out[i] , result[i] ) ;
- assert ( fabs ( out[i] - result[i] ) / fabs( out[i]) < 1e-6 ) ;
- }
-
-}
-
-static int testDiva (void) {
-
- printf("\n>>>> Float real Tests\n");
- /*sdivaTest();*/
- printf("\n>>>> Double real Tests\n");
- ddivaTest();
- return 0;
-}
-
-
-int main(void) {
- assert(testDiva() == 0);
- return 0;
-}
diff --git a/src/matrixOperations/division/testMatrixLDivision.c b/src/matrixOperations/division/testMatrixLDivision.c
index ac5d57bf..0dc9a6c6 100644
--- a/src/matrixOperations/division/testMatrixLDivision.c
+++ b/src/matrixOperations/division/testMatrixLDivision.c
@@ -10,18 +10,22 @@
*
*/
-#include "matrixDivision.h"
+
#include <assert.h>
#include <stdio.h>
+#include "matrixDivision.h"
+
+
+#define LINES 2
+#define COLUMNS1 2
+#define COLUMNS2 2
+
-#define LINES 20
-#define COLUMNS1 17
-#define COLUMNS2 12
static void dldivaTest ( void )
{
int i = 0 ;
-
- double in1[] =
+/* here the matrixes are linearized in the wrong way so need to transpose them */
+/* double in1[] =
{0.84155184263363481,0.26385784195736051,0.52570608118548989,0.54653349192813039,
0.62128817522898316,0.98085420625284314,0.74896081397309899,0.01432593585923314,
0.23678414756432176,0.70614896761253476,0.27255953708663583,0.06706167431548238,
@@ -239,11 +243,17 @@ static void dldivaTest ( void )
0.52872556238983059,0.57098625620273025,1.06212082783082407,0.59423966432136910,
0.34731906030365728,-1.35636228486699051,0.60027136476023302,-0.01921912593105047,
-0.26023872041264290,0.60546432820947238,-0.64774098495099597,1.07329929083653908}
-;
+; */
+
+ double in1[] = { 4 , 3 , 8 , 9 } ;
+ double in2[] = { 1 , 3 , 2 , 4 } ;
+ double result[] = { -1.25 , 0.75 , -1.166666666666666 ,0.83333333333333333 };
+
+
double out [(COLUMNS2)*(LINES)] ;
- dldiva ( in1 , LINES , COLUMNS1 , in2 , LINES , COLUMNS2 , out ) ;
+ dldiva( in1 , LINES , COLUMNS1 , in2 , LINES , COLUMNS2 , out ) ;
for ( i = 0 ; i < LINES *COLUMNS2 ; i++ )
{
printf ( "\n %d out : %e result : %e assert : %e \n" , i , out[i] , result[i] , fabs ( out[i] - result[i] ) / fabs( out[i]) ) ;
@@ -254,17 +264,22 @@ static void dldivaTest ( void )
}
+
+
+
static int testLDiva (void) {
- printf("\n>>>> Float real Tests\n");
+ printf("\n>>>> Left Tests\n");
/*sdivaTest();*/
- printf("\n>>>> Double real Tests\n");
+ printf("\n\t>>>> Double real Tests\n");
dldivaTest();
+
+
return 0;
}
-
int main(void) {
- assert(testLDiva() == 0);
+ assert(testLDiva () == 0);
return 0;
}
+
diff --git a/src/matrixOperations/division/testMatrixRDivision.c b/src/matrixOperations/division/testMatrixRDivision.c
index 3be29088..bc211664 100644
--- a/src/matrixOperations/division/testMatrixRDivision.c
+++ b/src/matrixOperations/division/testMatrixRDivision.c
@@ -10,13 +10,17 @@
*
*/
-#include "matrixDivision.h"
+
#include <assert.h>
#include <stdio.h>
+#include "matrixDivision.h"
#define LINES1 2
-#define LINES2 2
+#define LINES2 1
#define COLUMNS 2
+
+
+
/*
static void sdivaTest ( void )
{
@@ -218,40 +222,46 @@ float out[(COLUMNS)*(LINES2)] ;
*/
-/*
-static void ddivaTest ( void )
+
+
+static void drdivaTest ( void )
{
int i = 0 ;
- double in1[] = { 4 , 8 , 3 , 9 } ;
- double in2[] = { 1 , 2 , 3 , 4 } ;
- double result[] = { 4 , 0 , 7.5 ,-1.5 };
+/* double in1[] = { 4 , 3 , 8 , 9 } ;
+ double in2[] = { 1 , 3 , 2 , 4 } ;
+ double result[] = { 4 , 7.5 , 0 ,-1.5 };*/
+
+
+ double in1[] = { 1 ,3 ,2 ,4 } ;
+ double in2[] = { 1 , 2 } ;
+ double result[] = { 1 , 2.2 };
double out [(COLUMNS)*(LINES2)] ;
- drdiva ( in1 , LINES1 , COLUMNS , in2 , LINES2 , COLUMNS , out ) ;
+ drdiva ( in1 , LINES1 , COLUMNS , in2 , LINES2 , COLUMNS , out) ;
printf("\n\n\t>>>>>debut assert\n");
for ( i = 0 ; i < LINES2 *COLUMNS ; i++ )
{
printf ( "\n %d out : %e result : %e assert : %e \n" , i , out[i] , result[i] , fabs ( out[i] - result[i] ) / fabs( out[i]) ) ;
- assert ( fabs ( out[i] - result[i] ) / fabs( out[i]) < 100000000 ) ;
+ assert ( fabs ( out[i] - result[i] ) / fabs( out[i]) < 1e-15 ) ;
}
}
-*/
-static int testDiva (void) {
- printf("\n>>>> Float real Tests\n");
+static int testRDiva (void) {
+
+ printf("\n>>>> Right Tests\n");
/*sdivaTest();*/
- printf("\n>>>> Double real Tests\n");
- ddivaTest();
+ printf("\n\t>>>> Double real Tests\n");
+ drdivaTest();
return 0;
}
int main(void) {
- assert(testDiva() == 0);
+ assert(testRDiva () == 0);
return 0;
}
diff --git a/src/matrixOperations/division/zrdiva.c b/src/matrixOperations/division/zrdiva.c
index c1b56179..b5bd9dea 100644
--- a/src/matrixOperations/division/zrdiva.c
+++ b/src/matrixOperations/division/zrdiva.c
@@ -11,36 +11,132 @@
*/
/**** WARNING NOT WORK AT ALL FOR THE MOMENT ***/
-/*
-void zrdiva ( doubleComplex * in1, int it1,
- doubleComplex * in2, int it2,
- doubleComplex * out, int itOut, int size){
- int iIndex = 0; //Main loop index
- int iIndex1 = 0; //Loop index on left operand
- int iIndex2 = 0; //Loop index on right operand
- int iIndexOut = 0; //Lopp index on result matrix
+/**/
+
+#include "matrixDivision.h"
+#include "lapack.h"
+#include <stdio.h>
+#include <string.h>
+int zrdiva( doubleComplex* in1, int lines1, int columns1 ,
+ doubleComplex* in2, int lines2, int columns2 ,
+ doubleComplex* out )
+{
+
+ char cNorm = 0;
+ int iExit = 0;
+
+ /*temporary variables*/
+ int iWork = 0;
+ int iInfo = 0;
+ int iMax = 0;
+ double dblRcond = 0;
+
+ double dblEps = 0;
+ double dblAnorm = 0;
+
+
+ doubleComplex *poAf = NULL;
+ doubleComplex *poAt = NULL;
+ doubleComplex *poBt = NULL;
+ doubleComplex *poDwork = NULL;
+
+ int *pRank = NULL;
+ int *pIpiv = NULL;
+ int *pJpvt = NULL;
+ double *pRwork = NULL;
+
+ iWork = Max(2*columns2, Min(lines2, columns2) + Max(2 * Min(lines2, columns2), Max(lines2 + 1, Min(lines2, columns2) + lines1)));
+
+
+
+ /* Array allocations*/
+
+
+ poAf = (doubleComplex*)malloc(sizeof(doubleComplex) *(unsigned int) lines2 *(unsigned int) columns2);
+ poAt = (doubleComplex*)malloc(sizeof(doubleComplex) * (unsigned int)lines2 *(unsigned int) columns2);
+ poBt = (doubleComplex*)malloc(sizeof(doubleComplex) *(unsigned int) Max(lines2, columns2) *(unsigned int) lines1);
+
+ pRank = (int*)malloc(sizeof(int));
+ pIpiv = (int*)malloc(sizeof(int) *(unsigned int) columns2);
+ pJpvt = (int*)malloc(sizeof(int) *(unsigned int) lines2);
+ pRwork = (double*)malloc(sizeof(double) * 2 *(unsigned int) lines2);
+
+
+
+ cNorm = '1';
+ poDwork = (doubleComplex*)malloc(sizeof(doubleComplex) *(unsigned int) iWork);
+ dblEps = getRelativeMachinePrecision() ;
+ dblAnorm = C2F(zlange)(&cNorm, &lines2, &columns2, in2, &lines2, poDwork);
+ /*tranpose A and B*/
- for(iIndex = 0 ; iIndex < size ; iIndex++)
+ ztransposea(in2, lines2, columns2, poAt);
+ ztransposea(in1, lines1, columns2, poBt);
+
+ if(lines2 == columns2)
{
- out[iIndexOut] = zdivides( in1[iIndex1], in2[iIndex2] ) ;
-
- iIndexOut += itOut;
- iIndex1 += it1;
- iIndex2 += it2;
- }
-
-
-}
-
+ cNorm = 'F';
+ C2F(zlacpy)(&cNorm, &columns2, &columns2, poAt, &columns1, poAf, &columns2);
+ C2F(zgetrf)(&columns2, &columns2, poAf, &columns2, pIpiv, &iInfo);
+ if(iInfo == 0)
+ {
+ cNorm = '1';
+ C2F(zgecon)(&cNorm, &columns2, poAf, &columns2, &dblAnorm,
+ &dblRcond, poDwork, pRwork, &iInfo);
+ if(dblRcond > sqrt(dblEps))
+ {
+ cNorm = 'N';
+ C2F(zgetrs)(&cNorm, &columns2, &lines1, poAf, &columns2, pIpiv, poBt, &columns2, &iInfo);
+ ztransposea(poBt, columns2, lines2, out);
-void zldiva ( doubleComplex * in1, int it1,
- doubleComplex * in2, int it2,
- doubleComplex * out, int itOut, int size)) {
-f
-
- zrdiva ( in2 , it2 , in1 , it1 , out , itout , size );
+ iExit = 1;
+ }
+ }
+
+ }
+
+ if(iExit == 0)
+ {
+ dblRcond = sqrt(dblEps);
+ cNorm = 'F';
+ iMax = Max(lines2, columns2);
+ memset(pJpvt, 0x00,(unsigned int) sizeof(int) *(unsigned int) lines2);
+ C2F(zgelsy)(&columns2, &lines2, &lines1, poAt, &columns2, poBt, &iMax,
+ pJpvt, &dblRcond, &pRank[0], poDwork, &iWork, pRwork, &iInfo);
+
+ if(iInfo == 0)
+ {
+
+ /*// TransposeRealMatrix(pBt, lines1, lines2, _pdblRealOut, Max(lines1,columns1), lines2);
+
+ //Mega caca de la mort qui tue des ours a mains nues
+ //mais je ne sais pas comment le rendre "beau" :(*/
+ {
+ int i,j,ij,ji;
+ for(j = 0 ; j < lines2 ; j++)
+ {
+ for(i = 0 ; i < lines1 ; i++)
+ {
+ ij = i + j * lines1;
+ ji = j + i * Max(lines2, columns2);
+ out[ij] = DoubleComplex ( zreals( poBt[ji]) , -zimags ( poBt[ji]));
+ }
+ }
+ }
+ }
+ }
+
+
+
+ free(poAf);
+ free(poAt);
+ free(poBt);
+ free(pRank);
+ free(pIpiv);
+ free(pJpvt);
+ free(pRwork);
+ free(poDwork);
+ return 0;
}
-*/
generated by cgit (git 2.34.1) at 2025-03-21 16:04:08 +0000