Skip to content
Projects
Groups
Snippets
Help
Loading...
Help
Contribute to GitLab
Sign in
Toggle navigation
M
MarlinKimbra
Project
Project
Details
Activity
Cycle Analytics
Repository
Repository
Files
Commits
Branches
Tags
Contributors
Graph
Compare
Charts
Issues
0
Issues
0
List
Board
Labels
Milestones
Merge Requests
0
Merge Requests
0
CI / CD
CI / CD
Pipelines
Jobs
Schedules
Charts
Members
Members
Collapse sidebar
Close sidebar
Activity
Graph
Charts
Create a new issue
Jobs
Commits
Issue Boards
Open sidebar
machinery
MarlinKimbra
Commits
3653c5b0
Commit
3653c5b0
authored
Dec 04, 2015
by
MagoKimbra
Browse files
Options
Browse Files
Download
Email Patches
Plain Diff
Same fix
parent
078e021e
Changes
6
Show whitespace changes
Inline
Side-by-side
Showing
6 changed files
with
404 additions
and
735 deletions
+404
-735
Marlin_main.cpp
MarlinKimbra/Marlin_main.cpp
+3
-3
nextion_lcd.cpp
MarlinKimbra/nextion_lcd.cpp
+1
-1
nextion_lcd.h
MarlinKimbra/nextion_lcd.h
+1
-1
pins.h
MarlinKimbra/pins.h
+12
-5
qr_solve.cpp
MarlinKimbra/qr_solve.cpp
+369
-705
qr_solve.h
MarlinKimbra/qr_solve.h
+18
-20
No files found.
MarlinKimbra/Marlin_main.cpp
View file @
3653c5b0
...
@@ -353,7 +353,7 @@ unsigned long printer_usage_seconds;
...
@@ -353,7 +353,7 @@ unsigned long printer_usage_seconds;
#endif // FWRETRACT
#endif // FWRETRACT
#if
ENABLED(ULTIPANEL
) && HAS(POWER_SWITCH)
#if
(ENABLED(ULTIPANEL) || ENABLED(NEXTION)
) && HAS(POWER_SWITCH)
bool
powersupply
=
bool
powersupply
=
#if ENABLED(PS_DEFAULT_OFF)
#if ENABLED(PS_DEFAULT_OFF)
false
false
...
@@ -3692,13 +3692,13 @@ inline void gcode_G28() {
...
@@ -3692,13 +3692,13 @@ inline void gcode_G28() {
// raise extruder
// raise extruder
float
measured_z
,
float
measured_z
,
z_before
=
probePointCounter
?
Z_RAISE_BETWEEN_PROBINGS
+
current_position
[
Z_AXIS
]
:
Z_RAISE_BEFORE_PROBING
;
z_before
=
probePointCounter
?
Z_RAISE_BETWEEN_PROBINGS
+
current_position
[
Z_AXIS
]
:
Z_RAISE_BEFORE_PROBING
+
current_position
[
Z_AXIS
]
;
if
(
debugLevel
&
DEBUG_INFO
)
{
if
(
debugLevel
&
DEBUG_INFO
)
{
if
(
probePointCounter
)
if
(
probePointCounter
)
ECHO_LMV
(
DB
,
"z_before = (between) "
,
(
float
)(
Z_RAISE_BETWEEN_PROBINGS
+
current_position
[
Z_AXIS
]));
ECHO_LMV
(
DB
,
"z_before = (between) "
,
(
float
)(
Z_RAISE_BETWEEN_PROBINGS
+
current_position
[
Z_AXIS
]));
else
else
ECHO_LMV
(
DB
,
"z_before = (before) "
,
(
float
)
Z_RAISE_BEFORE_PROBING
);
ECHO_LMV
(
DB
,
"z_before = (before) "
,
(
float
)
Z_RAISE_BEFORE_PROBING
+
current_position
[
Z_AXIS
]
);
}
}
ProbeAction
act
;
ProbeAction
act
;
...
...
MarlinKimbra/nextion_lcd.cpp
View file @
3653c5b0
...
@@ -21,7 +21,7 @@
...
@@ -21,7 +21,7 @@
bool
PageInfo
=
false
;
bool
PageInfo
=
false
;
char
buffer
[
100
]
=
{
0
};
char
buffer
[
100
]
=
{
0
};
uint32_t
slidermaxval
=
20
;
uint32_t
slidermaxval
=
20
;
char
lcd_status_message
[
30
]
=
WELCOME_MSG
;
// worst case is kana with up to 3*LCD_WIDTH+1
char
lcd_status_message
[
30
]
=
WELCOME_MSG
;
uint8_t
lcd_status_message_level
=
0
;
uint8_t
lcd_status_message_level
=
0
;
static
millis_t
next_lcd_update_ms
;
static
millis_t
next_lcd_update_ms
;
...
...
MarlinKimbra/nextion_lcd.h
View file @
3653c5b0
...
@@ -4,7 +4,7 @@
...
@@ -4,7 +4,7 @@
#if ENABLED(NEXTION)
#if ENABLED(NEXTION)
#define LCD_UPDATE_INTERVAL 5000
#define LCD_UPDATE_INTERVAL 5000
void
Exit
Up
PopCallback
(
void
*
ptr
);
void
ExitPopCallback
(
void
*
ptr
);
void
setpagePopCallback
(
void
*
ptr
);
void
setpagePopCallback
(
void
*
ptr
);
void
hotPopCallback
(
void
*
ptr
);
void
hotPopCallback
(
void
*
ptr
);
void
sethotPopCallback
(
void
*
ptr
);
void
sethotPopCallback
(
void
*
ptr
);
...
...
MarlinKimbra/pins.h
View file @
3653c5b0
...
@@ -2297,15 +2297,22 @@
...
@@ -2297,15 +2297,22 @@
#define ORIG_E2_DIR_PIN 53
#define ORIG_E2_DIR_PIN 53
#define ORIG_E2_ENABLE_PIN 49
#define ORIG_E2_ENABLE_PIN 49
#define ORIG_E3_STEP_PIN 35
#define ORIG_E3_DIR_PIN 33
#define ORIG_E3_ENABLE_PIN 37
#define ORIG_E4_STEP_PIN 29
#define ORIG_E4_DIR_PIN 27
#define ORIG_E4_ENABLE_PIN 31
#define SDPOWER -1
#define SDPOWER -1
#define SDSS 4
#define SDSS 4
#define LED_PIN -1
#define LED_PIN -1
#define BEEPER_PIN 41
#define BEEPER_PIN 41
#define ORIG_FAN_PIN -1
#define ORIG_FAN_PIN 9
#define ORIG_FAN2_PIN 8
//#define CONTROLLERORIG_FAN_PIN 8 //Pin used for the fan to cool controller
#define PS_ON_PIN 40
#define PS_ON_PIN 40
...
...
MarlinKimbra/qr_solve.cpp
View file @
3653c5b0
#include "qr_solve.h"
#include "base.h"
#if ENABLED(AUTO_BED_LEVELING_FEATURE) && ENABLED(AUTO_BED_LEVELING_GRID)
#if ENABLED(AUTO_BED_LEVELING_FEATURE) && ENABLED(AUTO_BED_LEVELING_GRID)
#include "qr_solve.h"
#include <stdlib.h>
#include <math.h>
//# include "r8lib.h"
//# include "r8lib.h"
int
i4_min
(
int
i1
,
int
i2
)
int
i4_min
(
int
i1
,
int
i2
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -32,20 +34,10 @@ int i4_min ( int i1, int i2 )
...
@@ -32,20 +34,10 @@ int i4_min ( int i1, int i2 )
Output, int I4_MIN, the smaller of I1 and I2.
Output, int I4_MIN, the smaller of I1 and I2.
*/
*/
{
{
int
value
;
return
(
i1
<
i2
)
?
i1
:
i2
;
if
(
i1
<
i2
)
{
value
=
i1
;
}
else
{
value
=
i2
;
}
return
value
;
}
}
double
r8_epsilon
(
void
)
double
r8_epsilon
(
void
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -79,11 +71,10 @@ double r8_epsilon ( void )
...
@@ -79,11 +71,10 @@ double r8_epsilon ( void )
*/
*/
{
{
const
double
value
=
2.220446049250313E-016
;
const
double
value
=
2.220446049250313E-016
;
return
value
;
return
value
;
}
}
double
r8_max
(
double
x
,
double
y
)
double
r8_max
(
double
x
,
double
y
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -110,20 +101,10 @@ double r8_max ( double x, double y )
...
@@ -110,20 +101,10 @@ double r8_max ( double x, double y )
Output, double R8_MAX, the maximum of X and Y.
Output, double R8_MAX, the maximum of X and Y.
*/
*/
{
{
double
value
;
return
(
y
<
x
)
?
x
:
y
;
if
(
y
<
x
)
{
value
=
x
;
}
else
{
value
=
y
;
}
return
value
;
}
}
double
r8_abs
(
double
x
)
double
r8_abs
(
double
x
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -150,20 +131,10 @@ double r8_abs ( double x )
...
@@ -150,20 +131,10 @@ double r8_abs ( double x )
Output, double R8_ABS, the absolute value of X.
Output, double R8_ABS, the absolute value of X.
*/
*/
{
{
double
value
;
return
(
x
<
0.0
)
?
-
x
:
x
;
if
(
0.0
<=
x
)
{
value
=
+
x
;
}
else
{
value
=
-
x
;
}
return
value
;
}
}
double
r8_sign
(
double
x
)
double
r8_sign
(
double
x
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -190,20 +161,10 @@ double r8_sign ( double x )
...
@@ -190,20 +161,10 @@ double r8_sign ( double x )
Output, double R8_SIGN, the sign of X.
Output, double R8_SIGN, the sign of X.
*/
*/
{
{
double
value
;
return
(
x
<
0.0
)
?
-
1.0
:
1.0
;
if
(
x
<
0.0
)
{
value
=
-
1.0
;
}
else
{
value
=
+
1.0
;
}
return
value
;
}
}
double
r8mat_amax
(
int
m
,
int
n
,
double
a
[]
)
double
r8mat_amax
(
int
m
,
int
n
,
double
a
[]
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -239,26 +200,17 @@ double r8mat_amax ( int m, int n, double a[] )
...
@@ -239,26 +200,17 @@ double r8mat_amax ( int m, int n, double a[] )
Output, double R8MAT_AMAX, the maximum absolute value entry of A.
Output, double R8MAT_AMAX, the maximum absolute value entry of A.
*/
*/
{
{
int
i
;
double
value
=
r8_abs
(
a
[
0
+
0
*
m
]);
int
j
;
for
(
int
j
=
0
;
j
<
n
;
j
++
)
{
double
value
;
for
(
int
i
=
0
;
i
<
m
;
i
++
)
{
if
(
value
<
r8_abs
(
a
[
i
+
j
*
m
]))
value
=
r8_abs
(
a
[
0
+
0
*
m
]
);
value
=
r8_abs
(
a
[
i
+
j
*
m
]);
for
(
j
=
0
;
j
<
n
;
j
++
)
{
for
(
i
=
0
;
i
<
m
;
i
++
)
{
if
(
value
<
r8_abs
(
a
[
i
+
j
*
m
]
)
)
{
value
=
r8_abs
(
a
[
i
+
j
*
m
]
);
}
}
}
}
}
return
value
;
return
value
;
}
}
void
r8mat_copy
(
double
a2
[],
int
m
,
int
n
,
double
a1
[]
)
void
r8mat_copy
(
double
a2
[],
int
m
,
int
n
,
double
a1
[]
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -292,21 +244,15 @@ void r8mat_copy( double a2[], int m, int n, double a1[] )
...
@@ -292,21 +244,15 @@ void r8mat_copy( double a2[], int m, int n, double a1[] )
Output, double R8MAT_COPY_NEW[M*N], the copy of A1.
Output, double R8MAT_COPY_NEW[M*N], the copy of A1.
*/
*/
{
{
int
i
;
for
(
int
j
=
0
;
j
<
n
;
j
++
)
{
int
j
;
for
(
int
i
=
0
;
i
<
m
;
i
++
)
a2
[
i
+
j
*
m
]
=
a1
[
i
+
j
*
m
];
for
(
j
=
0
;
j
<
n
;
j
++
)
{
for
(
i
=
0
;
i
<
m
;
i
++
)
{
a2
[
i
+
j
*
m
]
=
a1
[
i
+
j
*
m
];
}
}
}
}
}
/******************************************************************************/
/******************************************************************************/
void
daxpy
(
int
n
,
double
da
,
double
dx
[],
int
incx
,
double
dy
[],
int
incy
)
void
daxpy
(
int
n
,
double
da
,
double
dx
[],
int
incx
,
double
dy
[],
int
incy
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -358,76 +304,46 @@ void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy )
...
@@ -358,76 +304,46 @@ void daxpy ( int n, double da, double dx[], int incx, double dy[], int incy )
Input, int INCY, the increment between successive entries of DY.
Input, int INCY, the increment between successive entries of DY.
*/
*/
{
{
int
i
;
if
(
n
<=
0
||
da
==
0.0
)
return
;
int
ix
;
int
iy
;
int
m
;
if
(
n
<=
0
)
int
i
,
ix
,
iy
,
m
;
{
/*
return
;
}
if
(
da
==
0.0
)
{
return
;
}
/*
Code for unequal increments or equal increments
Code for unequal increments or equal increments
not equal to 1.
not equal to 1.
*/
*/
if
(
incx
!=
1
||
incy
!=
1
)
if
(
incx
!=
1
||
incy
!=
1
)
{
{
if
(
0
<=
incx
)
if
(
0
<=
incx
)
{
ix
=
0
;
ix
=
0
;
}
else
else
{
ix
=
(
-
n
+
1
)
*
incx
;
ix
=
(
-
n
+
1
)
*
incx
;
if
(
0
<=
incy
)
}
if
(
0
<=
incy
)
{
iy
=
0
;
iy
=
0
;
}
else
else
{
iy
=
(
-
n
+
1
)
*
incy
;
iy
=
(
-
n
+
1
)
*
incy
;
for
(
i
=
0
;
i
<
n
;
i
++
)
{
}
for
(
i
=
0
;
i
<
n
;
i
++
)
{
dy
[
iy
]
=
dy
[
iy
]
+
da
*
dx
[
ix
];
dy
[
iy
]
=
dy
[
iy
]
+
da
*
dx
[
ix
];
ix
=
ix
+
incx
;
ix
=
ix
+
incx
;
iy
=
iy
+
incy
;
iy
=
iy
+
incy
;
}
}
}
}
/*
/*
Code for both increments equal to 1.
Code for both increments equal to 1.
*/
*/
else
else
{
{
m
=
n
%
4
;
m
=
n
%
4
;
for
(
i
=
0
;
i
<
m
;
i
++
)
for
(
i
=
0
;
i
<
m
;
i
++
)
{
dy
[
i
]
=
dy
[
i
]
+
da
*
dx
[
i
];
dy
[
i
]
=
dy
[
i
]
+
da
*
dx
[
i
];
}
for
(
i
=
m
;
i
<
n
;
i
=
i
+
4
)
{
for
(
i
=
m
;
i
<
n
;
i
=
i
+
4
)
{
dy
[
i
]
=
dy
[
i
]
+
da
*
dx
[
i
];
dy
[
i
]
=
dy
[
i
]
+
da
*
dx
[
i
];
dy
[
i
+
1
]
=
dy
[
i
+
1
]
+
da
*
dx
[
i
+
1
];
dy
[
i
+
1
]
=
dy
[
i
+
1
]
+
da
*
dx
[
i
+
1
];
dy
[
i
+
2
]
=
dy
[
i
+
2
]
+
da
*
dx
[
i
+
2
];
dy
[
i
+
2
]
=
dy
[
i
+
2
]
+
da
*
dx
[
i
+
2
];
dy
[
i
+
3
]
=
dy
[
i
+
3
]
+
da
*
dx
[
i
+
3
];
dy
[
i
+
3
]
=
dy
[
i
+
3
]
+
da
*
dx
[
i
+
3
];
}
}
}
}
return
;
}
}
/******************************************************************************/
/******************************************************************************/
double
ddot
(
int
n
,
double
dx
[],
int
incx
,
double
dy
[],
int
incy
)
double
ddot
(
int
n
,
double
dx
[],
int
incx
,
double
dy
[],
int
incy
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -479,75 +395,45 @@ double ddot ( int n, double dx[], int incx, double dy[], int incy )
...
@@ -479,75 +395,45 @@ double ddot ( int n, double dx[], int incx, double dy[], int incy )
entries of DX and DY.
entries of DX and DY.
*/
*/
{
{
double
dtemp
;
int
i
;
int
ix
;
int
iy
;
int
m
;
dtemp
=
0.0
;
if
(
n
<=
0
)
return
0.0
;
if
(
n
<=
0
)
int
i
,
m
;
{
double
dtemp
=
0.0
;
return
dtemp
;
}
/*
/*
Code for unequal increments or equal increments
Code for unequal increments or equal increments
not equal to 1.
not equal to 1.
*/
*/
if
(
incx
!=
1
||
incy
!=
1
)
if
(
incx
!=
1
||
incy
!=
1
)
{
{
int
ix
=
(
incx
>=
0
)
?
0
:
(
-
n
+
1
)
*
incx
,
if
(
0
<=
incx
)
iy
=
(
incy
>=
0
)
?
0
:
(
-
n
+
1
)
*
incy
;
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
ix
=
0
;
dtemp
+=
dx
[
ix
]
*
dy
[
iy
];
}
else
{
ix
=
(
-
n
+
1
)
*
incx
;
}
if
(
0
<=
incy
)
{
iy
=
0
;
}
else
{
iy
=
(
-
n
+
1
)
*
incy
;
}
for
(
i
=
0
;
i
<
n
;
i
++
)
{
dtemp
=
dtemp
+
dx
[
ix
]
*
dy
[
iy
];
ix
=
ix
+
incx
;
ix
=
ix
+
incx
;
iy
=
iy
+
incy
;
iy
=
iy
+
incy
;
}
}
}
}
/*
/*
Code for both increments equal to 1.
Code for both increments equal to 1.
*/
*/
else
else
{
{
m
=
n
%
5
;
m
=
n
%
5
;
for
(
i
=
0
;
i
<
m
;
i
++
)
for
(
i
=
0
;
i
<
m
;
i
++
)
dtemp
+=
dx
[
i
]
*
dy
[
i
];
{
for
(
i
=
m
;
i
<
n
;
i
=
i
+
5
)
{
dtemp
=
dtemp
+
dx
[
i
]
*
dy
[
i
];
dtemp
+=
dx
[
i
]
*
dy
[
i
]
}
+
dx
[
i
+
1
]
*
dy
[
i
+
1
]
+
dx
[
i
+
2
]
*
dy
[
i
+
2
]
for
(
i
=
m
;
i
<
n
;
i
=
i
+
5
)
+
dx
[
i
+
3
]
*
dy
[
i
+
3
]
{
+
dx
[
i
+
4
]
*
dy
[
i
+
4
];
dtemp
=
dtemp
+
dx
[
i
]
*
dy
[
i
]
+
dx
[
i
+
1
]
*
dy
[
i
+
1
]
+
dx
[
i
+
2
]
*
dy
[
i
+
2
]
+
dx
[
i
+
3
]
*
dy
[
i
+
3
]
+
dx
[
i
+
4
]
*
dy
[
i
+
4
];
}
}
}
}
return
dtemp
;
return
dtemp
;
}
}
/******************************************************************************/
/******************************************************************************/
double
dnrm2
(
int
n
,
double
x
[],
int
incx
)
double
dnrm2
(
int
n
,
double
x
[],
int
incx
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -594,54 +480,34 @@ double dnrm2 ( int n, double x[], int incx )
...
@@ -594,54 +480,34 @@ double dnrm2 ( int n, double x[], int incx )
Output, double DNRM2, the Euclidean norm of X.
Output, double DNRM2, the Euclidean norm of X.
*/
*/
{
{
double
absxi
;
int
i
;
int
ix
;
double
norm
;
double
norm
;
double
scale
;
if
(
n
<
1
||
incx
<
1
)
double
ssq
;
if
(
n
<
1
||
incx
<
1
)
{
norm
=
0.0
;
norm
=
0.0
;
}
else
if
(
n
==
1
)
else
if
(
n
==
1
)
norm
=
r8_abs
(
x
[
0
]);
{
else
{
norm
=
r8_abs
(
x
[
0
]
);
double
scale
=
0.0
,
ssq
=
1.0
;
}
int
ix
=
0
;
else
for
(
int
i
=
0
;
i
<
n
;
i
++
)
{
{
if
(
x
[
ix
]
!=
0.0
)
{
scale
=
0.0
;
double
absxi
=
r8_abs
(
x
[
ix
]);
ssq
=
1.0
;
if
(
scale
<
absxi
)
{
ix
=
0
;
ssq
=
1.0
+
ssq
*
(
scale
/
absxi
)
*
(
scale
/
absxi
);
for
(
i
=
0
;
i
<
n
;
i
++
)
{
if
(
x
[
ix
]
!=
0.0
)
{
absxi
=
r8_abs
(
x
[
ix
]
);
if
(
scale
<
absxi
)
{
ssq
=
1.0
+
ssq
*
(
scale
/
absxi
)
*
(
scale
/
absxi
);
scale
=
absxi
;
scale
=
absxi
;
}
}
else
else
{
ssq
=
ssq
+
(
absxi
/
scale
)
*
(
absxi
/
scale
);
ssq
=
ssq
+
(
absxi
/
scale
)
*
(
absxi
/
scale
);
}
}
ix
+=
incx
;
}
}
ix
=
ix
+
incx
;
norm
=
scale
*
sqrt
(
ssq
);
}
norm
=
scale
*
sqrt
(
ssq
);
}
}
return
norm
;
return
norm
;
}
}
/******************************************************************************/
/******************************************************************************/
void
dqrank
(
double
a
[],
int
lda
,
int
m
,
int
n
,
double
tol
,
int
*
kr
,
void
dqrank
(
double
a
[],
int
lda
,
int
m
,
int
n
,
double
tol
,
int
*
kr
,
int
jpvt
[],
double
qraux
[]
)
int
jpvt
[],
double
qraux
[]
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -715,39 +581,27 @@ void dqrank ( double a[], int lda, int m, int n, double tol, int *kr,
...
@@ -715,39 +581,27 @@ void dqrank ( double a[], int lda, int m, int n, double tol, int *kr,
the QR factorization.
the QR factorization.
*/
*/
{
{
int
i
;
int
j
;
int
job
;
int
k
;
double
work
[
n
];
double
work
[
n
];
for
(
i
=
0
;
i
<
n
;
i
++
)
for
(
int
i
=
0
;
i
<
n
;
i
++
)
{
jpvt
[
i
]
=
0
;
jpvt
[
i
]
=
0
;
}
job
=
1
;
int
job
=
1
;
dqrdc
(
a
,
lda
,
m
,
n
,
qraux
,
jpvt
,
work
,
job
);
dqrdc
(
a
,
lda
,
m
,
n
,
qraux
,
jpvt
,
work
,
job
);
*
kr
=
0
;
*
kr
=
0
;
k
=
i4_min
(
m
,
n
);
int
k
=
i4_min
(
m
,
n
);
for
(
int
j
=
0
;
j
<
k
;
j
++
)
{
for
(
j
=
0
;
j
<
k
;
j
++
)
if
(
r8_abs
(
a
[
j
+
j
*
lda
])
<=
tol
*
r8_abs
(
a
[
0
+
0
*
lda
]))
{
if
(
r8_abs
(
a
[
j
+
j
*
lda
]
)
<=
tol
*
r8_abs
(
a
[
0
+
0
*
lda
]
)
)
{
return
;
return
;
}
*
kr
=
j
+
1
;
*
kr
=
j
+
1
;
}
}
return
;
}
}
/******************************************************************************/
/******************************************************************************/
void
dqrdc
(
double
a
[],
int
lda
,
int
n
,
int
p
,
double
qraux
[],
int
jpvt
[],
void
dqrdc
(
double
a
[],
int
lda
,
int
n
,
int
p
,
double
qraux
[],
int
jpvt
[],
double
work
[],
int
job
)
double
work
[],
int
job
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -827,176 +681,121 @@ void dqrdc ( double a[], int lda, int n, int p, double qraux[], int jpvt[],
...
@@ -827,176 +681,121 @@ void dqrdc ( double a[], int lda, int n, int p, double qraux[], int jpvt[],
nonzero, pivoting is done.
nonzero, pivoting is done.
*/
*/
{
{
int
j
;
int
jp
;
int
jp
;
int
l
;
int
j
;
int
lup
;
int
lup
;
int
maxj
;
int
maxj
;
double
maxnrm
;
double
maxnrm
,
nrmxl
,
t
,
tt
;
double
nrmxl
;
int
pl
;
int
pu
;
int
swapj
;
double
t
;
double
tt
;
pl
=
1
;
int
pl
=
1
,
pu
=
0
;
pu
=
0
;
/*
/*
If pivoting is requested, rearrange the columns.
If pivoting is requested, rearrange the columns.
*/
*/
if
(
job
!=
0
)
if
(
job
!=
0
)
{
{
for
(
j
=
1
;
j
<=
p
;
j
++
)
{
for
(
j
=
1
;
j
<=
p
;
j
++
)
int
swapj
=
(
0
<
jpvt
[
j
-
1
]);
{
jpvt
[
j
-
1
]
=
(
jpvt
[
j
-
1
]
<
0
)
?
-
j
:
j
;
swapj
=
(
0
<
jpvt
[
j
-
1
]
);
if
(
swapj
)
{
if
(
j
!=
pl
)
if
(
jpvt
[
j
-
1
]
<
0
)
dswap
(
n
,
a
+
0
+
(
pl
-
1
)
*
lda
,
1
,
a
+
0
+
(
j
-
1
),
1
);
{
jpvt
[
j
-
1
]
=
jpvt
[
pl
-
1
];
jpvt
[
j
-
1
]
=
-
j
;
jpvt
[
pl
-
1
]
=
j
;
}
pl
++
;
else
{
jpvt
[
j
-
1
]
=
j
;
}
if
(
swapj
)
{
if
(
j
!=
pl
)
{
dswap
(
n
,
a
+
0
+
(
pl
-
1
)
*
lda
,
1
,
a
+
0
+
(
j
-
1
),
1
);
}
jpvt
[
j
-
1
]
=
jpvt
[
pl
-
1
];
jpvt
[
pl
-
1
]
=
j
;
pl
=
pl
+
1
;
}
}
}
}
pu
=
p
;
pu
=
p
;
for
(
j
=
p
;
1
<=
j
;
j
--
)
{
for
(
j
=
p
;
1
<=
j
;
j
--
)
if
(
jpvt
[
j
-
1
]
<
0
)
{
{
jpvt
[
j
-
1
]
=
-
jpvt
[
j
-
1
];
if
(
jpvt
[
j
-
1
]
<
0
)
if
(
j
!=
pu
)
{
{
dswap
(
n
,
a
+
0
+
(
pu
-
1
)
*
lda
,
1
,
a
+
0
+
(
j
-
1
)
*
lda
,
1
);
jpvt
[
j
-
1
]
=
-
jpvt
[
j
-
1
];
jp
=
jpvt
[
pu
-
1
];
jpvt
[
pu
-
1
]
=
jpvt
[
j
-
1
];
if
(
j
!=
pu
)
jpvt
[
j
-
1
]
=
jp
;
{
dswap
(
n
,
a
+
0
+
(
pu
-
1
)
*
lda
,
1
,
a
+
0
+
(
j
-
1
)
*
lda
,
1
);
jp
=
jpvt
[
pu
-
1
];
jpvt
[
pu
-
1
]
=
jpvt
[
j
-
1
];
jpvt
[
j
-
1
]
=
jp
;
}
}
pu
=
pu
-
1
;
pu
=
pu
-
1
;
}
}
}
}
}
}
/*
/*
Compute the norms of the free columns.
Compute the norms of the free columns.
*/
*/
for
(
j
=
pl
;
j
<=
pu
;
j
++
)
for
(
j
=
pl
;
j
<=
pu
;
j
++
)
{
qraux
[
j
-
1
]
=
dnrm2
(
n
,
a
+
0
+
(
j
-
1
)
*
lda
,
1
);
qraux
[
j
-
1
]
=
dnrm2
(
n
,
a
+
0
+
(
j
-
1
)
*
lda
,
1
);
for
(
j
=
pl
;
j
<=
pu
;
j
++
)
}
work
[
j
-
1
]
=
qraux
[
j
-
1
];
/*
for
(
j
=
pl
;
j
<=
pu
;
j
++
)
{
work
[
j
-
1
]
=
qraux
[
j
-
1
];
}
/*
Perform the Householder reduction of A.
Perform the Householder reduction of A.
*/
*/
lup
=
i4_min
(
n
,
p
);
lup
=
i4_min
(
n
,
p
);
for
(
int
l
=
1
;
l
<=
lup
;
l
++
)
{
for
(
l
=
1
;
l
<=
lup
;
l
++
)
/*
{
/*
Bring the column of largest norm into the pivot position.
Bring the column of largest norm into the pivot position.
*/
*/
if
(
pl
<=
l
&&
l
<
pu
)
if
(
pl
<=
l
&&
l
<
pu
)
{
{
maxnrm
=
0.0
;
maxnrm
=
0.0
;
maxj
=
l
;
maxj
=
l
;
for
(
j
=
l
;
j
<=
pu
;
j
++
)
for
(
j
=
l
;
j
<=
pu
;
j
++
)
{
{
if
(
maxnrm
<
qraux
[
j
-
1
])
{
if
(
maxnrm
<
qraux
[
j
-
1
]
)
maxnrm
=
qraux
[
j
-
1
];
{
maxnrm
=
qraux
[
j
-
1
];
maxj
=
j
;
maxj
=
j
;
}
}
}
}
if
(
maxj
!=
l
)
{
if
(
maxj
!=
l
)
dswap
(
n
,
a
+
0
+
(
l
-
1
)
*
lda
,
1
,
a
+
0
+
(
maxj
-
1
)
*
lda
,
1
);
{
qraux
[
maxj
-
1
]
=
qraux
[
l
-
1
];
dswap
(
n
,
a
+
0
+
(
l
-
1
)
*
lda
,
1
,
a
+
0
+
(
maxj
-
1
)
*
lda
,
1
);
work
[
maxj
-
1
]
=
work
[
l
-
1
];
qraux
[
maxj
-
1
]
=
qraux
[
l
-
1
];
jp
=
jpvt
[
maxj
-
1
];
work
[
maxj
-
1
]
=
work
[
l
-
1
];
jpvt
[
maxj
-
1
]
=
jpvt
[
l
-
1
];
jp
=
jpvt
[
maxj
-
1
];
jpvt
[
l
-
1
]
=
jp
;
jpvt
[
maxj
-
1
]
=
jpvt
[
l
-
1
];
jpvt
[
l
-
1
]
=
jp
;
}
}
}
}
/*
/*
Compute the Householder transformation for column L.
Compute the Householder transformation for column L.
*/
*/
qraux
[
l
-
1
]
=
0.0
;
qraux
[
l
-
1
]
=
0.0
;
if
(
l
!=
n
)
{
if
(
l
!=
n
)
nrmxl
=
dnrm2
(
n
-
l
+
1
,
a
+
l
-
1
+
(
l
-
1
)
*
lda
,
1
);
{
if
(
nrmxl
!=
0.0
)
{
nrmxl
=
dnrm2
(
n
-
l
+
1
,
a
+
l
-
1
+
(
l
-
1
)
*
lda
,
1
);
if
(
a
[
l
-
1
+
(
l
-
1
)
*
lda
]
!=
0.0
)
nrmxl
=
nrmxl
*
r8_sign
(
a
[
l
-
1
+
(
l
-
1
)
*
lda
]);
if
(
nrmxl
!=
0.0
)
dscal
(
n
-
l
+
1
,
1.0
/
nrmxl
,
a
+
l
-
1
+
(
l
-
1
)
*
lda
,
1
);
{
a
[
l
-
1
+
(
l
-
1
)
*
lda
]
=
1.0
+
a
[
l
-
1
+
(
l
-
1
)
*
lda
];
if
(
a
[
l
-
1
+
(
l
-
1
)
*
lda
]
!=
0.0
)
/*
{
nrmxl
=
nrmxl
*
r8_sign
(
a
[
l
-
1
+
(
l
-
1
)
*
lda
]
);
}
dscal
(
n
-
l
+
1
,
1.0
/
nrmxl
,
a
+
l
-
1
+
(
l
-
1
)
*
lda
,
1
);
a
[
l
-
1
+
(
l
-
1
)
*
lda
]
=
1.0
+
a
[
l
-
1
+
(
l
-
1
)
*
lda
];
/*
Apply the transformation to the remaining columns, updating the norms.
Apply the transformation to the remaining columns, updating the norms.
*/
*/
for
(
j
=
l
+
1
;
j
<=
p
;
j
++
)
for
(
j
=
l
+
1
;
j
<=
p
;
j
++
)
{
{
t
=
-
ddot
(
n
-
l
+
1
,
a
+
l
-
1
+
(
l
-
1
)
*
lda
,
1
,
a
+
l
-
1
+
(
j
-
1
)
*
lda
,
1
)
t
=
-
ddot
(
n
-
l
+
1
,
a
+
l
-
1
+
(
l
-
1
)
*
lda
,
1
,
a
+
l
-
1
+
(
j
-
1
)
*
lda
,
1
)
/
a
[
l
-
1
+
(
l
-
1
)
*
lda
];
/
a
[
l
-
1
+
(
l
-
1
)
*
lda
];
daxpy
(
n
-
l
+
1
,
t
,
a
+
l
-
1
+
(
l
-
1
)
*
lda
,
1
,
a
+
l
-
1
+
(
j
-
1
)
*
lda
,
1
);
daxpy
(
n
-
l
+
1
,
t
,
a
+
l
-
1
+
(
l
-
1
)
*
lda
,
1
,
a
+
l
-
1
+
(
j
-
1
)
*
lda
,
1
);
if
(
pl
<=
j
&&
j
<=
pu
)
{
if
(
qraux
[
j
-
1
]
!=
0.0
)
{
if
(
pl
<=
j
&&
j
<=
pu
)
tt
=
1.0
-
pow
(
r8_abs
(
a
[
l
-
1
+
(
j
-
1
)
*
lda
])
/
qraux
[
j
-
1
],
2
);
{
tt
=
r8_max
(
tt
,
0.0
);
if
(
qraux
[
j
-
1
]
!=
0.0
)
{
tt
=
1.0
-
pow
(
r8_abs
(
a
[
l
-
1
+
(
j
-
1
)
*
lda
]
)
/
qraux
[
j
-
1
],
2
);
tt
=
r8_max
(
tt
,
0.0
);
t
=
tt
;
t
=
tt
;
tt
=
1.0
+
0.05
*
tt
*
pow
(
qraux
[
j
-
1
]
/
work
[
j
-
1
],
2
);
tt
=
1.0
+
0.05
*
tt
*
pow
(
qraux
[
j
-
1
]
/
work
[
j
-
1
],
2
);
if
(
tt
!=
1.0
)
if
(
tt
!=
1.0
)
qraux
[
j
-
1
]
=
qraux
[
j
-
1
]
*
sqrt
(
t
);
{
else
{
qraux
[
j
-
1
]
=
qraux
[
j
-
1
]
*
sqrt
(
t
);
qraux
[
j
-
1
]
=
dnrm2
(
n
-
l
,
a
+
l
+
(
j
-
1
)
*
lda
,
1
);
}
work
[
j
-
1
]
=
qraux
[
j
-
1
];
else
{
qraux
[
j
-
1
]
=
dnrm2
(
n
-
l
,
a
+
l
+
(
j
-
1
)
*
lda
,
1
);
work
[
j
-
1
]
=
qraux
[
j
-
1
];
}
}
}
}
}
}
}
}
/*
/*
Save the transformation.
Save the transformation.
*/
*/
qraux
[
l
-
1
]
=
a
[
l
-
1
+
(
l
-
1
)
*
lda
];
qraux
[
l
-
1
]
=
a
[
l
-
1
+
(
l
-
1
)
*
lda
];
a
[
l
-
1
+
(
l
-
1
)
*
lda
]
=
-
nrmxl
;
a
[
l
-
1
+
(
l
-
1
)
*
lda
]
=
-
nrmxl
;
}
}
}
}
}
}
return
;
}
}
/******************************************************************************/
/******************************************************************************/
int
dqrls
(
double
a
[],
int
lda
,
int
m
,
int
n
,
double
tol
,
int
*
kr
,
double
b
[],
int
dqrls
(
double
a
[],
int
lda
,
int
m
,
int
n
,
double
tol
,
int
*
kr
,
double
b
[],
double
x
[],
double
rsd
[],
int
jpvt
[],
double
qraux
[],
int
itask
)
double
x
[],
double
rsd
[],
int
jpvt
[],
double
qraux
[],
int
itask
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -1104,9 +903,7 @@ int dqrls ( double a[], int lda, int m, int n, double tol, int *kr, double b[],
...
@@ -1104,9 +903,7 @@ int dqrls ( double a[], int lda, int m, int n, double tol, int *kr, double b[],
*/
*/
{
{
int
ind
;
int
ind
;
if
(
lda
<
m
)
{
if
(
lda
<
m
)
{
/*fprintf ( stderr, "\n" );
/*fprintf ( stderr, "\n" );
fprintf ( stderr, "DQRLS - Fatal error!\n" );
fprintf ( stderr, "DQRLS - Fatal error!\n" );
fprintf ( stderr, " LDA < M.\n" );*/
fprintf ( stderr, " LDA < M.\n" );*/
...
@@ -1114,8 +911,7 @@ int dqrls ( double a[], int lda, int m, int n, double tol, int *kr, double b[],
...
@@ -1114,8 +911,7 @@ int dqrls ( double a[], int lda, int m, int n, double tol, int *kr, double b[],
return
ind
;
return
ind
;
}
}
if
(
n
<=
0
)
if
(
n
<=
0
)
{
{
/*fprintf ( stderr, "\n" );
/*fprintf ( stderr, "\n" );
fprintf ( stderr, "DQRLS - Fatal error!\n" );
fprintf ( stderr, "DQRLS - Fatal error!\n" );
fprintf ( stderr, " N <= 0.\n" );*/
fprintf ( stderr, " N <= 0.\n" );*/
...
@@ -1123,8 +919,7 @@ int dqrls ( double a[], int lda, int m, int n, double tol, int *kr, double b[],
...
@@ -1123,8 +919,7 @@ int dqrls ( double a[], int lda, int m, int n, double tol, int *kr, double b[],
return
ind
;
return
ind
;
}
}
if
(
itask
<
1
)
if
(
itask
<
1
)
{
{
/*fprintf ( stderr, "\n" );
/*fprintf ( stderr, "\n" );
fprintf ( stderr, "DQRLS - Fatal error!\n" );
fprintf ( stderr, "DQRLS - Fatal error!\n" );
fprintf ( stderr, " ITASK < 1.\n" );*/
fprintf ( stderr, " ITASK < 1.\n" );*/
...
@@ -1133,24 +928,21 @@ int dqrls ( double a[], int lda, int m, int n, double tol, int *kr, double b[],
...
@@ -1133,24 +928,21 @@ int dqrls ( double a[], int lda, int m, int n, double tol, int *kr, double b[],
}
}
ind
=
0
;
ind
=
0
;
/*
/*
Factor the matrix.
Factor the matrix.
*/
*/
if
(
itask
==
1
)
if
(
itask
==
1
)
{
dqrank
(
a
,
lda
,
m
,
n
,
tol
,
kr
,
jpvt
,
qraux
);
dqrank
(
a
,
lda
,
m
,
n
,
tol
,
kr
,
jpvt
,
qraux
);
/*
}
/*
Solve the least-squares problem.
Solve the least-squares problem.
*/
*/
dqrlss
(
a
,
lda
,
m
,
n
,
*
kr
,
b
,
x
,
rsd
,
jpvt
,
qraux
);
dqrlss
(
a
,
lda
,
m
,
n
,
*
kr
,
b
,
x
,
rsd
,
jpvt
,
qraux
);
return
ind
;
return
ind
;
}
}
/******************************************************************************/
/******************************************************************************/
void
dqrlss
(
double
a
[],
int
lda
,
int
m
,
int
n
,
int
kr
,
double
b
[],
double
x
[],
void
dqrlss
(
double
a
[],
int
lda
,
int
m
,
int
n
,
int
kr
,
double
b
[],
double
x
[],
double
rsd
[],
int
jpvt
[],
double
qraux
[]
)
double
rsd
[],
int
jpvt
[],
double
qraux
[]
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -1230,45 +1022,37 @@ void dqrlss ( double a[], int lda, int m, int n, int kr, double b[], double x[],
...
@@ -1230,45 +1022,37 @@ void dqrlss ( double a[], int lda, int m, int n, int kr, double b[], double x[],
int
k
;
int
k
;
double
t
;
double
t
;
if
(
kr
!=
0
)
if
(
kr
!=
0
)
{
{
job
=
110
;
job
=
110
;
info
=
dqrsl
(
a
,
lda
,
m
,
kr
,
qraux
,
b
,
rsd
,
rsd
,
x
,
rsd
,
rsd
,
job
);
info
=
dqrsl
(
a
,
lda
,
m
,
kr
,
qraux
,
b
,
rsd
,
rsd
,
x
,
rsd
,
rsd
,
job
);
UNUSED
(
info
);
}
}
for
(
i
=
0
;
i
<
n
;
i
++
)
for
(
i
=
0
;
i
<
n
;
i
++
)
{
jpvt
[
i
]
=
-
jpvt
[
i
];
jpvt
[
i
]
=
-
jpvt
[
i
];
}
for
(
i
=
kr
;
i
<
n
;
i
++
)
for
(
i
=
kr
;
i
<
n
;
i
++
)
{
x
[
i
]
=
0.0
;
x
[
i
]
=
0.0
;
}
for
(
j
=
1
;
j
<=
n
;
j
++
)
for
(
j
=
1
;
j
<=
n
;
j
++
)
{
{
if
(
jpvt
[
j
-
1
]
<=
0
)
{
if
(
jpvt
[
j
-
1
]
<=
0
)
k
=
-
jpvt
[
j
-
1
];
{
jpvt
[
j
-
1
]
=
k
;
k
=
-
jpvt
[
j
-
1
];
jpvt
[
j
-
1
]
=
k
;
while
(
k
!=
j
)
while
(
k
!=
j
)
{
{
t
=
x
[
j
-
1
];
t
=
x
[
j
-
1
];
x
[
j
-
1
]
=
x
[
k
-
1
];
x
[
j
-
1
]
=
x
[
k
-
1
];
x
[
k
-
1
]
=
t
;
x
[
k
-
1
]
=
t
;
jpvt
[
k
-
1
]
=
-
jpvt
[
k
-
1
];
jpvt
[
k
-
1
]
=
-
jpvt
[
k
-
1
];
k
=
jpvt
[
k
-
1
];
k
=
jpvt
[
k
-
1
];
}
}
}
}
}
}
return
;
}
}
/******************************************************************************/
/******************************************************************************/
int
dqrsl
(
double
a
[],
int
lda
,
int
n
,
int
k
,
double
qraux
[],
double
y
[],
int
dqrsl
(
double
a
[],
int
lda
,
int
n
,
int
k
,
double
qraux
[],
double
y
[],
double
qy
[],
double
qty
[],
double
b
[],
double
rsd
[],
double
ab
[],
int
job
)
double
qy
[],
double
qty
[],
double
b
[],
double
rsd
[],
double
ab
[],
int
job
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -1418,217 +1202,151 @@ int dqrsl ( double a[], int lda, int n, int k, double qraux[], double y[],
...
@@ -1418,217 +1202,151 @@ int dqrsl ( double a[], int lda, int n, int k, double qraux[], double y[],
int
ju
;
int
ju
;
double
t
;
double
t
;
double
temp
;
double
temp
;
/*
/*
Set INFO flag.
Set INFO flag.
*/
*/
info
=
0
;
info
=
0
;
/*
Determine what is to be computed.
*/
cqy
=
(
job
/
10000
!=
0
);
cqty
=
(
(
job
%
10000
)
!=
0
);
cb
=
(
(
job
%
1000
)
/
100
!=
0
);
cr
=
(
(
job
%
100
)
/
10
!=
0
);
cab
=
(
(
job
%
10
)
!=
0
);
ju
=
i4_min
(
k
,
n
-
1
);
/*
/*
Determine what is to be computed.
*/
cqy
=
(
job
/
10000
!=
0
);
cqty
=
((
job
%
10000
)
!=
0
);
cb
=
((
job
%
1000
)
/
100
!=
0
);
cr
=
((
job
%
100
)
/
10
!=
0
);
cab
=
((
job
%
10
)
!=
0
);
ju
=
i4_min
(
k
,
n
-
1
);
/*
Special action when N = 1.
Special action when N = 1.
*/
*/
if
(
ju
==
0
)
if
(
ju
==
0
)
{
{
if
(
cqy
)
if
(
cqy
)
{
qy
[
0
]
=
y
[
0
];
qy
[
0
]
=
y
[
0
];
}
if
(
cqty
)
if
(
cqty
)
{
qty
[
0
]
=
y
[
0
];
qty
[
0
]
=
y
[
0
];
}
if
(
cab
)
if
(
cab
)
{
ab
[
0
]
=
y
[
0
];
ab
[
0
]
=
y
[
0
];
}
if
(
cb
)
{
if
(
a
[
0
+
0
*
lda
]
==
0.0
)
if
(
cb
)
{
if
(
a
[
0
+
0
*
lda
]
==
0.0
)
{
info
=
1
;
info
=
1
;
}
else
else
{
b
[
0
]
=
y
[
0
]
/
a
[
0
+
0
*
lda
];
b
[
0
]
=
y
[
0
]
/
a
[
0
+
0
*
lda
];
}
}
}
if
(
cr
)
if
(
cr
)
{
rsd
[
0
]
=
0.0
;
rsd
[
0
]
=
0.0
;
}
return
info
;
return
info
;
}
}
/*
/*
Set up to compute QY or QTY.
Set up to compute QY or QTY.
*/
*/
if
(
cqy
)
if
(
cqy
)
{
{
for
(
i
=
1
;
i
<=
n
;
i
++
)
for
(
i
=
1
;
i
<=
n
;
i
++
)
qy
[
i
-
1
]
=
y
[
i
-
1
];
{
qy
[
i
-
1
]
=
y
[
i
-
1
];
}
}
if
(
cqty
)
{
for
(
i
=
1
;
i
<=
n
;
i
++
)
qty
[
i
-
1
]
=
y
[
i
-
1
];
}
}
/*
if
(
cqty
)
{
for
(
i
=
1
;
i
<=
n
;
i
++
)
{
qty
[
i
-
1
]
=
y
[
i
-
1
];
}
}
/*
Compute QY.
Compute QY.
*/
*/
if
(
cqy
)
if
(
cqy
)
{
{
for
(
jj
=
1
;
jj
<=
ju
;
jj
++
)
{
for
(
jj
=
1
;
jj
<=
ju
;
jj
++
)
{
j
=
ju
-
jj
+
1
;
j
=
ju
-
jj
+
1
;
if
(
qraux
[
j
-
1
]
!=
0.0
)
{
if
(
qraux
[
j
-
1
]
!=
0.0
)
temp
=
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
{
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
qraux
[
j
-
1
];
temp
=
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
t
=
-
ddot
(
n
-
j
+
1
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
qy
+
j
-
1
,
1
)
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
qraux
[
j
-
1
];
daxpy
(
n
-
j
+
1
,
t
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
qy
+
j
-
1
,
1
);
t
=
-
ddot
(
n
-
j
+
1
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
qy
+
j
-
1
,
1
)
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
temp
;
daxpy
(
n
-
j
+
1
,
t
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
qy
+
j
-
1
,
1
);
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
temp
;
}
}
}
}
}
}
/*
/*
Compute Q'*Y.
Compute Q'*Y.
*/
*/
if
(
cqty
)
if
(
cqty
)
{
{
for
(
j
=
1
;
j
<=
ju
;
j
++
)
{
for
(
j
=
1
;
j
<=
ju
;
j
++
)
if
(
qraux
[
j
-
1
]
!=
0.0
)
{
{
temp
=
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
if
(
qraux
[
j
-
1
]
!=
0.0
)
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
qraux
[
j
-
1
];
{
t
=
-
ddot
(
n
-
j
+
1
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
qty
+
j
-
1
,
1
)
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
temp
=
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
daxpy
(
n
-
j
+
1
,
t
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
qty
+
j
-
1
,
1
);
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
qraux
[
j
-
1
];
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
temp
;
t
=
-
ddot
(
n
-
j
+
1
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
qty
+
j
-
1
,
1
)
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
daxpy
(
n
-
j
+
1
,
t
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
qty
+
j
-
1
,
1
);
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
temp
;
}
}
}
}
}
}
/*
/*
Set up to compute B, RSD, or AB.
Set up to compute B, RSD, or AB.
*/
*/
if
(
cb
)
if
(
cb
)
{
{
for
(
i
=
1
;
i
<=
k
;
i
++
)
for
(
i
=
1
;
i
<=
k
;
i
++
)
b
[
i
-
1
]
=
qty
[
i
-
1
];
{
b
[
i
-
1
]
=
qty
[
i
-
1
];
}
}
if
(
cab
)
{
for
(
i
=
1
;
i
<=
k
;
i
++
)
{
ab
[
i
-
1
]
=
qty
[
i
-
1
];
}
}
}
if
(
cab
)
{
if
(
cr
&&
k
<
n
)
for
(
i
=
1
;
i
<=
k
;
i
++
)
{
ab
[
i
-
1
]
=
qty
[
i
-
1
];
for
(
i
=
k
+
1
;
i
<=
n
;
i
++
)
{
rsd
[
i
-
1
]
=
qty
[
i
-
1
];
}
}
if
(
cr
&&
k
<
n
)
{
for
(
i
=
k
+
1
;
i
<=
n
;
i
++
)
rsd
[
i
-
1
]
=
qty
[
i
-
1
];
}
}
if
(
cab
&&
k
+
1
<=
n
)
{
if
(
cab
&&
k
+
1
<=
n
)
for
(
i
=
k
+
1
;
i
<=
n
;
i
++
)
{
ab
[
i
-
1
]
=
0.0
;
for
(
i
=
k
+
1
;
i
<=
n
;
i
++
)
{
ab
[
i
-
1
]
=
0.0
;
}
}
if
(
cr
)
{
for
(
i
=
1
;
i
<=
k
;
i
++
)
rsd
[
i
-
1
]
=
0.0
;
}
}
/*
if
(
cr
)
{
for
(
i
=
1
;
i
<=
k
;
i
++
)
{
rsd
[
i
-
1
]
=
0.0
;
}
}
/*
Compute B.
Compute B.
*/
*/
if
(
cb
)
if
(
cb
)
{
{
for
(
jj
=
1
;
jj
<=
k
;
jj
++
)
{
for
(
jj
=
1
;
jj
<=
k
;
jj
++
)
{
j
=
k
-
jj
+
1
;
j
=
k
-
jj
+
1
;
if
(
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
==
0.0
)
{
if
(
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
==
0.0
)
{
info
=
j
;
info
=
j
;
break
;
break
;
}
}
b
[
j
-
1
]
=
b
[
j
-
1
]
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
b
[
j
-
1
]
=
b
[
j
-
1
]
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
if
(
j
!=
1
)
{
t
=
-
b
[
j
-
1
];
if
(
j
!=
1
)
daxpy
(
j
-
1
,
t
,
a
+
0
+
(
j
-
1
)
*
lda
,
1
,
b
,
1
);
{
t
=
-
b
[
j
-
1
];
daxpy
(
j
-
1
,
t
,
a
+
0
+
(
j
-
1
)
*
lda
,
1
,
b
,
1
);
}
}
}
}
}
}
/*
/*
Compute RSD or AB as required.
Compute RSD or AB as required.
*/
*/
if
(
cr
||
cab
)
if
(
cr
||
cab
)
{
{
for
(
jj
=
1
;
jj
<=
ju
;
jj
++
)
{
for
(
jj
=
1
;
jj
<=
ju
;
jj
++
)
{
j
=
ju
-
jj
+
1
;
j
=
ju
-
jj
+
1
;
if
(
qraux
[
j
-
1
]
!=
0.0
)
{
if
(
qraux
[
j
-
1
]
!=
0.0
)
temp
=
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
{
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
qraux
[
j
-
1
];
temp
=
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
if
(
cr
)
{
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
qraux
[
j
-
1
];
t
=
-
ddot
(
n
-
j
+
1
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
rsd
+
j
-
1
,
1
)
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
if
(
cr
)
daxpy
(
n
-
j
+
1
,
t
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
rsd
+
j
-
1
,
1
);
{
t
=
-
ddot
(
n
-
j
+
1
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
rsd
+
j
-
1
,
1
)
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
daxpy
(
n
-
j
+
1
,
t
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
rsd
+
j
-
1
,
1
);
}
}
if
(
cab
)
{
if
(
cab
)
t
=
-
ddot
(
n
-
j
+
1
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
ab
+
j
-
1
,
1
)
{
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
t
=
-
ddot
(
n
-
j
+
1
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
ab
+
j
-
1
,
1
)
daxpy
(
n
-
j
+
1
,
t
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
ab
+
j
-
1
,
1
);
/
a
[
j
-
1
+
(
j
-
1
)
*
lda
];
daxpy
(
n
-
j
+
1
,
t
,
a
+
j
-
1
+
(
j
-
1
)
*
lda
,
1
,
ab
+
j
-
1
,
1
);
}
}
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
temp
;
a
[
j
-
1
+
(
j
-
1
)
*
lda
]
=
temp
;
}
}
}
}
}
}
return
info
;
return
info
;
}
}
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
void
dscal
(
int
n
,
double
sa
,
double
x
[],
int
incx
)
void
dscal
(
int
n
,
double
sa
,
double
x
[],
int
incx
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -1675,50 +1393,35 @@ void dscal ( int n, double sa, double x[], int incx )
...
@@ -1675,50 +1393,35 @@ void dscal ( int n, double sa, double x[], int incx )
int
ix
;
int
ix
;
int
m
;
int
m
;
if
(
n
<=
0
)
if
(
n
<=
0
)
return
;
{
}
else
if
(
incx
==
1
)
{
m
=
n
%
5
;
for
(
i
=
0
;
i
<
m
;
i
++
)
if
(
incx
==
1
)
{
{
m
=
n
%
5
;
for
(
i
=
0
;
i
<
m
;
i
++
)
x
[
i
]
=
sa
*
x
[
i
];
x
[
i
]
=
sa
*
x
[
i
];
}
for
(
i
=
m
;
i
<
n
;
i
=
i
+
5
)
{
for
(
i
=
m
;
i
<
n
;
i
=
i
+
5
)
{
x
[
i
]
=
sa
*
x
[
i
];
x
[
i
]
=
sa
*
x
[
i
];
x
[
i
+
1
]
=
sa
*
x
[
i
+
1
];
x
[
i
+
1
]
=
sa
*
x
[
i
+
1
];
x
[
i
+
2
]
=
sa
*
x
[
i
+
2
];
x
[
i
+
2
]
=
sa
*
x
[
i
+
2
];
x
[
i
+
3
]
=
sa
*
x
[
i
+
3
];
x
[
i
+
3
]
=
sa
*
x
[
i
+
3
];
x
[
i
+
4
]
=
sa
*
x
[
i
+
4
];
x
[
i
+
4
]
=
sa
*
x
[
i
+
4
];
}
}
}
}
else
else
{
{
if
(
0
<=
incx
)
if
(
0
<=
incx
)
{
ix
=
0
;
ix
=
0
;
}
else
else
{
ix
=
(
-
n
+
1
)
*
incx
;
ix
=
(
-
n
+
1
)
*
incx
;
for
(
i
=
0
;
i
<
n
;
i
++
)
{
}
for
(
i
=
0
;
i
<
n
;
i
++
)
{
x
[
ix
]
=
sa
*
x
[
ix
];
x
[
ix
]
=
sa
*
x
[
ix
];
ix
=
ix
+
incx
;
ix
=
ix
+
incx
;
}
}
}
}
return
;
}
}
/******************************************************************************/
/******************************************************************************/
void
dswap
(
int
n
,
double
x
[],
int
incx
,
double
y
[],
int
incy
)
void
dswap
(
int
n
,
double
x
[],
int
incx
,
double
y
[],
int
incy
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -1763,79 +1466,47 @@ void dswap ( int n, double x[], int incx, double y[], int incy )
...
@@ -1763,79 +1466,47 @@ void dswap ( int n, double x[], int incx, double y[], int incy )
Input, int INCY, the increment between successive elements of Y.
Input, int INCY, the increment between successive elements of Y.
*/
*/
{
{
int
i
;
if
(
n
<=
0
)
return
;
int
ix
;
int
iy
;
int
i
,
ix
,
iy
,
m
;
int
m
;
double
temp
;
double
temp
;
if
(
n
<=
0
)
if
(
incx
==
1
&&
incy
==
1
)
{
{
}
else
if
(
incx
==
1
&&
incy
==
1
)
{
m
=
n
%
3
;
m
=
n
%
3
;
for
(
i
=
0
;
i
<
m
;
i
++
)
{
for
(
i
=
0
;
i
<
m
;
i
++
)
{
temp
=
x
[
i
];
temp
=
x
[
i
];
x
[
i
]
=
y
[
i
];
x
[
i
]
=
y
[
i
];
y
[
i
]
=
temp
;
y
[
i
]
=
temp
;
}
}
for
(
i
=
m
;
i
<
n
;
i
=
i
+
3
)
{
for
(
i
=
m
;
i
<
n
;
i
=
i
+
3
)
{
temp
=
x
[
i
];
temp
=
x
[
i
];
x
[
i
]
=
y
[
i
];
x
[
i
]
=
y
[
i
];
y
[
i
]
=
temp
;
y
[
i
]
=
temp
;
temp
=
x
[
i
+
1
];
temp
=
x
[
i
+
1
];
x
[
i
+
1
]
=
y
[
i
+
1
];
x
[
i
+
1
]
=
y
[
i
+
1
];
y
[
i
+
1
]
=
temp
;
y
[
i
+
1
]
=
temp
;
temp
=
x
[
i
+
2
];
x
[
i
+
2
]
=
y
[
i
+
2
];
temp
=
x
[
i
+
2
];
y
[
i
+
2
]
=
temp
;
x
[
i
+
2
]
=
y
[
i
+
2
];
}
y
[
i
+
2
]
=
temp
;
}
}
else
{
}
ix
=
(
incx
>=
0
)
?
0
:
(
-
n
+
1
)
*
incx
;
else
iy
=
(
incy
>=
0
)
?
0
:
(
-
n
+
1
)
*
incy
;
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
if
(
0
<=
incx
)
{
ix
=
0
;
}
else
{
ix
=
(
-
n
+
1
)
*
incx
;
}
if
(
0
<=
incy
)
{
iy
=
0
;
}
else
{
iy
=
(
-
n
+
1
)
*
incy
;
}
for
(
i
=
0
;
i
<
n
;
i
++
)
{
temp
=
x
[
ix
];
temp
=
x
[
ix
];
x
[
ix
]
=
y
[
iy
];
x
[
ix
]
=
y
[
iy
];
y
[
iy
]
=
temp
;
y
[
iy
]
=
temp
;
ix
=
ix
+
incx
;
ix
=
ix
+
incx
;
iy
=
iy
+
incy
;
iy
=
iy
+
incy
;
}
}
}
}
return
;
}
}
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
/******************************************************************************/
void
qr_solve
(
double
x
[],
int
m
,
int
n
,
double
a
[],
double
b
[]
)
void
qr_solve
(
double
x
[],
int
m
,
int
n
,
double
a
[],
double
b
[]
)
/******************************************************************************/
/******************************************************************************/
/*
/*
...
@@ -1885,22 +1556,15 @@ void qr_solve ( double x[], int m, int n, double a[], double b[] )
...
@@ -1885,22 +1556,15 @@ void qr_solve ( double x[], int m, int n, double a[], double b[] )
Output, double QR_SOLVE[N], the least squares solution.
Output, double QR_SOLVE[N], the least squares solution.
*/
*/
{
{
double
a_qr
[
n
*
m
];
double
a_qr
[
n
*
m
],
qraux
[
n
],
r
[
m
],
tol
;
int
ind
;
int
ind
,
itask
,
jpvt
[
n
],
kr
,
lda
;
int
itask
;
int
jpvt
[
n
];
r8mat_copy
(
a_qr
,
m
,
n
,
a
);
int
kr
;
int
lda
;
double
qraux
[
n
];
double
r
[
m
];
double
tol
;
r8mat_copy
(
a_qr
,
m
,
n
,
a
);
lda
=
m
;
lda
=
m
;
tol
=
r8_epsilon
(
)
/
r8mat_amax
(
m
,
n
,
a_qr
);
tol
=
r8_epsilon
()
/
r8mat_amax
(
m
,
n
,
a_qr
);
itask
=
1
;
itask
=
1
;
ind
=
dqrls
(
a_qr
,
lda
,
m
,
n
,
tol
,
&
kr
,
b
,
x
,
r
,
jpvt
,
qraux
,
itask
);
ind
=
dqrls
(
a_qr
,
lda
,
m
,
n
,
tol
,
&
kr
,
b
,
x
,
r
,
jpvt
,
qraux
,
itask
);
UNUSED
(
ind
);
}
}
/******************************************************************************/
/******************************************************************************/
...
...
MarlinKimbra/qr_solve.h
View file @
3653c5b0
#i
f ENABLED(AUTO_BED_LEVELING_GRID)
#i
nclude "base.h"
#ifndef QR_SOLVE_H
#if ENABLED(AUTO_BED_LEVELING_GRID)
#define QR_SOLVE_H
void
daxpy
(
int
n
,
double
da
,
double
dx
[],
int
incx
,
double
dy
[],
int
incy
);
void
daxpy
(
int
n
,
double
da
,
double
dx
[],
int
incx
,
double
dy
[],
int
incy
);
double
ddot
(
int
n
,
double
dx
[],
int
incx
,
double
dy
[],
int
incy
);
double
ddot
(
int
n
,
double
dx
[],
int
incx
,
double
dy
[],
int
incy
);
double
dnrm2
(
int
n
,
double
x
[],
int
incx
);
double
dnrm2
(
int
n
,
double
x
[],
int
incx
);
void
dqrank
(
double
a
[],
int
lda
,
int
m
,
int
n
,
double
tol
,
int
*
kr
,
void
dqrank
(
double
a
[],
int
lda
,
int
m
,
int
n
,
double
tol
,
int
*
kr
,
int
jpvt
[],
double
qraux
[]
);
int
jpvt
[],
double
qraux
[]
);
void
dqrdc
(
double
a
[],
int
lda
,
int
n
,
int
p
,
double
qraux
[],
int
jpvt
[],
void
dqrdc
(
double
a
[],
int
lda
,
int
n
,
int
p
,
double
qraux
[],
int
jpvt
[],
double
work
[],
int
job
);
double
work
[],
int
job
);
int
dqrls
(
double
a
[],
int
lda
,
int
m
,
int
n
,
double
tol
,
int
*
kr
,
double
b
[],
int
dqrls
(
double
a
[],
int
lda
,
int
m
,
int
n
,
double
tol
,
int
*
kr
,
double
b
[],
double
x
[],
double
rsd
[],
int
jpvt
[],
double
qraux
[],
int
itask
);
double
x
[],
double
rsd
[],
int
jpvt
[],
double
qraux
[],
int
itask
);
void
dqrlss
(
double
a
[],
int
lda
,
int
m
,
int
n
,
int
kr
,
double
b
[],
double
x
[],
void
dqrlss
(
double
a
[],
int
lda
,
int
m
,
int
n
,
int
kr
,
double
b
[],
double
x
[],
double
rsd
[],
int
jpvt
[],
double
qraux
[]
);
double
rsd
[],
int
jpvt
[],
double
qraux
[]
);
int
dqrsl
(
double
a
[],
int
lda
,
int
n
,
int
k
,
double
qraux
[],
double
y
[],
int
dqrsl
(
double
a
[],
int
lda
,
int
n
,
int
k
,
double
qraux
[],
double
y
[],
double
qy
[],
double
qty
[],
double
b
[],
double
rsd
[],
double
ab
[],
int
job
);
double
qy
[],
double
qty
[],
double
b
[],
double
rsd
[],
double
ab
[],
int
job
);
void
dscal
(
int
n
,
double
sa
,
double
x
[],
int
incx
);
void
dscal
(
int
n
,
double
sa
,
double
x
[],
int
incx
);
void
dswap
(
int
n
,
double
x
[],
int
incx
,
double
y
[],
int
incy
);
void
dswap
(
int
n
,
double
x
[],
int
incx
,
double
y
[],
int
incy
);
void
qr_solve
(
double
x
[],
int
m
,
int
n
,
double
a
[],
double
b
[]
);
void
qr_solve
(
double
x
[],
int
m
,
int
n
,
double
a
[],
double
b
[]
);
#endif
#endif
#endif
\ No newline at end of file
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment