1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
|
/*
* 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"
#include "lapack.h"
void u16ldivma (uint16* in1, int lines1, int columns1 ,
uint16* in2, int lines2, int columns2 ,
uint16* out ){
char cNorm = 0;
int iExit = 0;
/*temporary variables*/
int iWork = 0;
int iInfo = 0;
int iMax = 0;
uint16 dblRcond = 0;
uint16 dblEps = 0;
uint16 dblAnorm = 0;
uint16 *pAf = NULL;
uint16 *pXb = NULL;
uint16 *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));
lines2 = 0 ;
/* Array allocations*/
pAf = (uint16*)malloc(sizeof(uint16) * (unsigned int) lines1 * (unsigned int) columns1);
pXb = (uint16*)malloc(sizeof(uint16) * (unsigned int) max(lines1,columns1) * (unsigned int) columns2);
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);
cNorm = '1';
pDwork = (uint16*)malloc(sizeof(uint16) *(unsigned int)iWork);
dblEps = getRelativeMachinePrecision() ;
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 = '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, &lines1, 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(iInfo == 0)
{
cNorm = 'F';
C2F(dlacpy)(&cNorm, &columns1, &columns2, pXb, &iMax, out, &columns1);
}
}
free(pAf);
free(pXb);
free(pRank);
free(pIpiv);
free(pJpvt);
free(pIwork);
free(pDwork);
}
|
