From db464f35f5a10b58d9ed1085e0b462689adee583 Mon Sep 17 00:00:00 2001
From: Siddhesh Wani
Date: Mon, 25 May 2015 14:46:31 +0530
Subject: Original Version
---
src/c/matrixOperations/division/zrdivma.c | 165 ++++++++++++++++++++++++++++++
1 file changed, 165 insertions(+)
create mode 100644 src/c/matrixOperations/division/zrdivma.c
(limited to 'src/c/matrixOperations/division/zrdivma.c')
diff --git a/src/c/matrixOperations/division/zrdivma.c b/src/c/matrixOperations/division/zrdivma.c
new file mode 100644
index 00000000..73e69f4c
--- /dev/null
+++ b/src/c/matrixOperations/division/zrdivma.c
@@ -0,0 +1,165 @@
+/*
+ * 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 <stdlib.h>
+#include <string.h>
+#include "matrixDivision.h"
+
+void zrdivma( 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);
+
+ /* hermitian tranpose A (transpose + conjugate )*/
+
+ ztransposea(in2, lines2, columns2, poAt);
+ zconja ( poAt , lines2*columns2 , poAt );
+ {
+ int i,j,ij,ji;
+ for(j = 0 ; j < lines1 ; j++)
+ {
+ for(i = 0 ; i < columns2 ; i++)
+ {
+ ij = i + j * max(lines2, columns2);
+ ji = j + i * lines1;
+
+ poBt[ij] = DoubleComplex (zreals ( in1[ji] ) , - zimags ( in1[ji] ) );
+ }
+ }
+ }
+
+ if(lines2 == columns2)
+ {
+ 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);
+ {
+ 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]));
+
+ }
+ }
+ }
+
+ 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);
+
+}
+
+
--
cgit