summaryrefslogtreecommitdiff
path: root/spectro/conv.c
diff options
context:
space:
mode:
authorJörg Frings-Fürst <debian@jff-webhosting.net>2016-10-02 19:24:58 +0200
committerJörg Frings-Fürst <debian@jff-webhosting.net>2016-10-02 19:24:58 +0200
commit3db384424bd7398ffbb7a355cab8f15f3add009f (patch)
tree4536961c62454aca3ac87ee88229e4d20c0d44fa /spectro/conv.c
parentd479dd1aab1c1cb907932c6595b0ef33523fc797 (diff)
New upstream version 1.9.1+repackupstream/1.9.1+repack
Diffstat (limited to 'spectro/conv.c')
-rw-r--r--spectro/conv.c659
1 files changed, 71 insertions, 588 deletions
diff --git a/spectro/conv.c b/spectro/conv.c
index 94023db..284a30c 100644
--- a/spectro/conv.c
+++ b/spectro/conv.c
@@ -64,7 +64,7 @@
#include "conv.h"
#include "icoms.h"
-#ifdef __APPLE__
+#ifdef UNIX_APPLE
//#include <stdbool.h>
#include <sys/sysctl.h>
#include <sys/param.h>
@@ -81,7 +81,7 @@
#if MAC_OS_X_VERSION_MIN_REQUIRED >= 1060
# include <objc/objc-auto.h>
#endif
-#endif /* __APPLE__ */
+#endif /* UNIX_APPLE */
#undef DEBUG
@@ -131,6 +131,8 @@ int next_con_char(void) {
return c;
}
+#ifdef NEVER // Don't need this now.
+
/* Horrible hack to poll stdin when we're not interactive. */
/* This has the drawback that the char and returm must be */
/* written in one operation for the character to be recognised - */
@@ -143,9 +145,12 @@ static int th_read_char(void *pp) {
DWORD bread;
if ((stdinh = GetStdHandle(STD_INPUT_HANDLE)) == INVALID_HANDLE_VALUE) {
+ *rp = 0; /* We've started */
return 0;
}
+ *rp = 0; /* We've started */
+
if (ReadFile(stdinh, buf, 1, &bread, NULL)
&& bread == 1
&& buf[0] != '\r' && buf[0] != '\n') {
@@ -154,32 +159,65 @@ static int th_read_char(void *pp) {
return 0;
}
+#endif /* NEVER */
/* If there is one, return the next character from the keyboard, else return 0 */
/* (If not_interactive, always returns 0) */
int poll_con_char(void) {
if (not_interactive) { /* Can't assume that it's the console */
+
+#ifdef NEVER // Use better approach below.
athread *getch_thread = NULL;
- char c = 0;
+ volatile char c = 0xff;
/* This is pretty horrible. The problem is that we can't use */
/* any of MSWin's async file read functions, because we */
/* have no way of ensuring that the STD_INPUT_HANDLE has been */
/* opened with FILE_FLAG_OVERLAPPED. Used a thread instead... */
- /* ReOpenFile() would fix this, but it's not available in WinXP, only Visa+ :-( */
- if ((getch_thread = new_athread(th_read_char, &c)) != NULL) {
+ /* ReOpenFile() would in theory fix this, but it's not available in WinXP, only Visa+, */
+ /* and aparently doesn't work on stdin anyway! :-( */
+ if ((getch_thread = new_athread(th_read_char, (char *)&c)) != NULL) {
HANDLE stdinh;
if ((stdinh = GetStdHandle(STD_INPUT_HANDLE)) == INVALID_HANDLE_VALUE) {
return 0;
}
- Sleep(100); /* We just hope 1 msec is enough for the thread to start */
- CancelIo(stdinh);
+ /* Wait for the thread to start */
+ while (c == 0xff) {
+ Sleep(10); /* We just hope 1 msec is enough for the thread to start */
+ }
+ Sleep(10); /* Give it time to read */
+ CancelIo(stdinh); /* May not work since ReadFile() is on a different thread ? */
getch_thread->del(getch_thread);
return c;
}
+#else /* ! NEVER */
+ /* This approach is very flakey from the console, but seems */
+ /* to work reliably when operated progromatically. */
+ HANDLE stdinh;
+ char buf[10] = { 0 }, c;
+ DWORD bread;
+
+ if ((stdinh = GetStdHandle(STD_INPUT_HANDLE)) == INVALID_HANDLE_VALUE) {
+ return 0;
+ }
+// printf("Waiting\n");
+ if (WaitForSingleObject(stdinh, 1) == WAIT_OBJECT_0) {
+// printf("stdin signalled\n");
+
+// FlushFileBuffers(stdinh);
+// FlushConsoleInputBuffer(stdinh);
+ if (ReadFile(stdinh, buf, 1, &bread, NULL)) {
+ int i;
+// fprintf(stderr,"Read %d chars 0x%x 0x%x 0x%x\n",bread,buf[0],buf[1], buf[2]);
+ if (buf[0] != '\r' && buf[0] != '\n')
+ return buf[0];
+ return 0;
+ }
+ }
+#endif /* !NEVER */
return 0;
}
@@ -192,10 +230,23 @@ int poll_con_char(void) {
}
/* Suck all characters from the keyboard */
-/* (If not_interactive, does nothing) */
+/* (If not_interactive, does nothing ?) */
void empty_con_chars(void) {
if (not_interactive) {
+ HANDLE stdinh;
+ char buf[100] = { 0 }, c;
+ DWORD bread;
+
+ if ((stdinh = GetStdHandle(STD_INPUT_HANDLE)) == INVALID_HANDLE_VALUE)
+ return;
+ for (;;) {
+ if (WaitForSingleObject(stdinh, 1) == WAIT_OBJECT_0) {
+ ReadFile(stdinh, buf, 100, &bread, NULL);
+ } else {
+ break;
+ }
+ }
return;
}
@@ -211,51 +262,6 @@ void sleep(unsigned int secs) {
Sleep(secs * 1000);
}
-/* Sleep for the given number of msec */
-void msec_sleep(unsigned int msec) {
- Sleep(msec);
-}
-
-/* Return the current time in msec since */
-/* the first invokation of msec_time() */
-/* (Is this based on timeGetTime() ? ) */
-unsigned int msec_time() {
- unsigned int rv;
- static unsigned int startup = 0;
-
- rv = GetTickCount();
- if (startup == 0)
- startup = rv;
-
- return rv - startup;
-}
-
-/* Return the current time in usec */
-/* since the first invokation of usec_time() */
-/* Return -1.0 if not available */
-double usec_time() {
- double rv;
- LARGE_INTEGER val;
- static double scale = 0.0;
- static LARGE_INTEGER startup;
-
- if (scale == 0.0) {
- if (QueryPerformanceFrequency(&val) == 0)
- return -1.0;
- scale = 1000000.0/val.QuadPart;
- QueryPerformanceCounter(&val);
- startup.QuadPart = val.QuadPart;
-
- } else {
- QueryPerformanceCounter(&val);
- }
- val.QuadPart -= startup.QuadPart;
-
- rv = val.QuadPart * scale;
-
- return rv;
-}
-
static athread *beep_thread = NULL;
static int beep_delay;
static int beep_freq;
@@ -552,126 +558,6 @@ void empty_con_chars(void) {
tcflush(STDIN_FILENO, TCIFLUSH);
}
-/* Sleep for the given number of msec */
-/* (Note that OS X 10.9+ App Nap can wreck this, unless */
-/* it is turned off.) */
-void msec_sleep(unsigned int msec) {
-#ifdef NEVER
- if (msec > 1000) {
- unsigned int secs;
- secs = msec / 1000;
- msec = msec % 1000;
- sleep(secs);
- }
- usleep(msec * 1000);
-#else
- struct timespec ts;
-
- ts.tv_sec = msec / 1000;
- ts.tv_nsec = (msec % 1000) * 1000000;
- nanosleep(&ts, NULL);
-#endif
-}
-
-
-#if defined(__APPLE__) && !defined(CLOCK_MONOTONIC)
-
-#include <mach/mach_time.h>
-
-unsigned int msec_time() {
- mach_timebase_info_data_t timebase;
- static uint64_t startup = 0;
- uint64_t time;
- double msec;
-
- time = mach_absolute_time();
- if (startup == 0)
- startup = time;
-
- mach_timebase_info(&timebase);
- time -= startup;
- msec = ((double)time * (double)timebase.numer)/((double)timebase.denom * 1e6);
-
- return (unsigned int)floor(msec + 0.5);
-}
-
-/* Return the current time in usec */
-/* since the first invokation of usec_time() */
-double usec_time() {
- mach_timebase_info_data_t timebase;
- static uint64_t startup = 0;
- uint64_t time;
- double usec;
-
- time = mach_absolute_time();
- if (startup == 0)
- startup = time;
-
- mach_timebase_info(&timebase);
- time -= startup;
- usec = ((double)time * (double)timebase.numer)/((double)timebase.denom * 1e3);
-
- return usec;
-}
-
-#else
-
-/* Return the current time in msec */
-/* since the first invokation of msec_time() */
-unsigned int msec_time() {
- unsigned int rv;
- static struct timespec startup = { 0, 0 };
- struct timespec cv;
-
- clock_gettime(CLOCK_MONOTONIC, &cv);
-
- /* Set time to 0 on first invocation */
- if (startup.tv_sec == 0 && startup.tv_nsec == 0)
- startup = cv;
-
- /* Subtract, taking care of carry */
- cv.tv_sec -= startup.tv_sec;
- if (startup.tv_nsec > cv.tv_nsec) {
- cv.tv_sec--;
- cv.tv_nsec += 1000000000;
- }
- cv.tv_nsec -= startup.tv_nsec;
-
- /* Convert nsec to msec */
- rv = cv.tv_sec * 1000 + cv.tv_nsec / 1000000;
-
- return rv;
-}
-
-/* Return the current time in usec */
-/* since the first invokation of usec_time() */
-double usec_time() {
- double rv;
- static struct timespec startup = { 0, 0 };
- struct timespec cv;
-
- clock_gettime(CLOCK_MONOTONIC, &cv);
-
- /* Set time to 0 on first invocation */
- if (startup.tv_sec == 0 && startup.tv_nsec == 0)
- startup = cv;
-
- /* Subtract, taking care of carry */
- cv.tv_sec -= startup.tv_sec;
- if (startup.tv_nsec > cv.tv_nsec) {
- cv.tv_sec--;
- cv.tv_nsec += 1000000000;
- }
- cv.tv_nsec -= startup.tv_nsec;
-
- /* Convert to usec */
- rv = cv.tv_sec * 1000000.0 + cv.tv_nsec/1000;
-
- return rv;
-}
-
-#endif
-
/* - - - - - - - - - - - - - - - - - - - - - - - - */
int acond_timedwait_imp(pthread_cond_t *cond, pthread_mutex_t *lock, int msec) {
@@ -699,7 +585,7 @@ int acond_timedwait_imp(pthread_cond_t *cond, pthread_mutex_t *lock, int msec) {
/* Set the current threads priority */
int set_interactive_priority() {
-#ifdef __APPLE__
+#ifdef UNIX_APPLE
#ifdef NEVER
int rv = 0;
struct task_category_policy tcatpolicy;
@@ -725,7 +611,7 @@ int set_interactive_priority() {
// a1logd(g_log, 8, "set_interactive_priority: set to important got %d\n",rv);
return rv;
#endif /* NEVER */
-#else /* !APPLE */
+#else /* !UNIX_APPLE */
int rv;
struct sched_param param;
param.sched_priority = 32;
@@ -734,11 +620,11 @@ int set_interactive_priority() {
rv = pthread_setschedparam(pthread_self(), SCHED_FIFO, &param);
// a1logd(g_log, 8, "set_interactive_priority: set got %d\n",rv);
return rv;
-#endif /* !APPLE */
+#endif /* !UNIX_APPLE */
}
int set_normal_priority() {
-#ifdef __APPLE__
+#ifdef UNIX_APPLE
#ifdef NEVER
int rv = 0;
struct task_category_policy tcatpolicy;
@@ -763,7 +649,7 @@ int set_normal_priority() {
// a1logd(g_log, 8, "set_normal_priority: set to standard got %d\n",rv);
return rv;
#endif /* NEVER */
-#else /* !APPLE */
+#else /* !UNIX_APPLE */
struct sched_param param;
param.sched_priority = 0;
int rv;
@@ -771,7 +657,7 @@ int set_normal_priority() {
rv = pthread_setschedparam(pthread_self(), SCHED_OTHER, &param);
// a1logd(g_log, 8, "set_normal_priority: reset got %d\n",rv);
return rv;
-#endif /* !APPLE */
+#endif /* !UNIX_APPLE */
}
#endif /* NEVER */
@@ -786,7 +672,7 @@ static int beep_msec;
static int delayed_beep(void *pp) {
msec_sleep(beep_delay);
a1logd(g_log,8, "msec_beep activate\n");
-#ifdef __APPLE__
+#ifdef UNIX_APPLE
# if __MAC_OS_X_VERSION_MAX_ALLOWED >= 1050
AudioServicesPlayAlertSound(kUserPreferredAlert);
# else
@@ -812,7 +698,7 @@ void msec_beep(int delay, int freq, int msec) {
a1logw(g_log, "msec_beep: Delayed beep failed to create thread\n");
} else {
a1logd(g_log,8, "msec_beep activate\n");
-#ifdef __APPLE__
+#ifdef UNIX_APPLE
# if __MAC_OS_X_VERSION_MAX_ALLOWED >= 1050
AudioServicesPlayAlertSound(kUserPreferredAlert);
# else
@@ -980,7 +866,7 @@ int create_parent_directories(char *path) {
#endif /* defined(UNIX) */
/* - - - - - - - - - - - - - - - - - - - - - - - - */
-#if defined(__APPLE__) || defined(NT)
+#if defined(UNIX_APPLE) || defined(NT)
/* Thread to monitor and kill the named processes */
static int th_kkill_nprocess(void *pp) {
@@ -1109,7 +995,7 @@ int kill_nprocess(char **pname, a1log *log) {
#endif /* NT */
-#if defined(__APPLE__)
+#if defined(UNIX_APPLE)
/* Kill a list of named process. */
/* Kill the first process found, then return */
@@ -1188,411 +1074,8 @@ int kill_nprocess(char **pname, a1log *log) {
return rv;
}
-#endif /* __APPLE__ */
-
-#endif /* __APPLE__ || NT */
-
-/* ============================================================= */
-
-/* A very small subset of icclib, copied to here. */
-/* This is just enough to support the standalone instruments */
-#ifdef SALONEINSTLIB
-
-sa_XYZNumber sa_D50 = {
- 0.9642, 1.0000, 0.8249
-};
-
-sa_XYZNumber sa_D65 = {
- 0.9505, 1.0000, 1.0890
-};
-
-void sa_SetUnity3x3(double mat[3][3]) {
- int i, j;
- for (j = 0; j < 3; j++) {
- for (i = 0; i < 3; i++) {
- if (i == j)
- mat[j][i] = 1.0;
- else
- mat[j][i] = 0.0;
- }
- }
-
-}
-
-void sa_Cpy3x3(double dst[3][3], double src[3][3]) {
- int i, j;
-
- for (j = 0; j < 3; j++) {
- for (i = 0; i < 3; i++)
- dst[j][i] = src[j][i];
- }
-}
-
-void sa_MulBy3x3(double out[3], double mat[3][3], double in[3]) {
- double tt[3];
-
- tt[0] = mat[0][0] * in[0] + mat[0][1] * in[1] + mat[0][2] * in[2];
- tt[1] = mat[1][0] * in[0] + mat[1][1] * in[1] + mat[1][2] * in[2];
- tt[2] = mat[2][0] * in[0] + mat[2][1] * in[1] + mat[2][2] * in[2];
+#endif /* UNIX_APPLE */
- out[0] = tt[0];
- out[1] = tt[1];
- out[2] = tt[2];
-}
-
-void sa_Mul3x3_2(double dst[3][3], double src1[3][3], double src2[3][3]) {
- int i, j, k;
- double td[3][3]; /* Temporary dest */
-
- for (j = 0; j < 3; j++) {
- for (i = 0; i < 3; i++) {
- double tt = 0.0;
- for (k = 0; k < 3; k++)
- tt += src1[j][k] * src2[k][i];
- td[j][i] = tt;
- }
- }
-
- /* Copy result out */
- for (j = 0; j < 3; j++)
- for (i = 0; i < 3; i++)
- dst[j][i] = td[j][i];
-}
-
-
-/* Matrix Inversion by Richard Carling from "Graphics Gems", Academic Press, 1990 */
-#define det2x2(a, b, c, d) (a * d - b * c)
-
-static void adjoint(
-double out[3][3],
-double in[3][3]
-) {
- double a1, a2, a3, b1, b2, b3, c1, c2, c3;
-
- /* assign to individual variable names to aid */
- /* selecting correct values */
-
- a1 = in[0][0]; b1 = in[0][1]; c1 = in[0][2];
- a2 = in[1][0]; b2 = in[1][1]; c2 = in[1][2];
- a3 = in[2][0]; b3 = in[2][1]; c3 = in[2][2];
-
- /* row column labeling reversed since we transpose rows & columns */
-
- out[0][0] = det2x2(b2, b3, c2, c3);
- out[1][0] = - det2x2(a2, a3, c2, c3);
- out[2][0] = det2x2(a2, a3, b2, b3);
-
- out[0][1] = - det2x2(b1, b3, c1, c3);
- out[1][1] = det2x2(a1, a3, c1, c3);
- out[2][1] = - det2x2(a1, a3, b1, b3);
-
- out[0][2] = det2x2(b1, b2, c1, c2);
- out[1][2] = - det2x2(a1, a2, c1, c2);
- out[2][2] = det2x2(a1, a2, b1, b2);
-}
-
-static double sa_Det3x3(double in[3][3]) {
- double a1, a2, a3, b1, b2, b3, c1, c2, c3;
- double ans;
-
- a1 = in[0][0]; b1 = in[0][1]; c1 = in[0][2];
- a2 = in[1][0]; b2 = in[1][1]; c2 = in[1][2];
- a3 = in[2][0]; b3 = in[2][1]; c3 = in[2][2];
-
- ans = a1 * det2x2(b2, b3, c2, c3)
- - b1 * det2x2(a2, a3, c2, c3)
- + c1 * det2x2(a2, a3, b2, b3);
- return ans;
-}
-
-#define SA__SMALL_NUMBER 1.e-8
-
-int sa_Inverse3x3(double out[3][3], double in[3][3]) {
- int i, j;
- double det;
-
- /* calculate the 3x3 determinant
- * if the determinant is zero,
- * then the inverse matrix is not unique.
- */
- det = sa_Det3x3(in);
-
- if ( fabs(det) < SA__SMALL_NUMBER)
- return 1;
-
- /* calculate the adjoint matrix */
- adjoint(out, in);
-
- /* scale the adjoint matrix to get the inverse */
- for (i = 0; i < 3; i++)
- for(j = 0; j < 3; j++)
- out[i][j] /= det;
- return 0;
-}
-
-#undef SA__SMALL_NUMBER
-#undef det2x2
-
-/* - - - - - - - - - - - - - - - - - - - - - - - - */
-/* Transpose a 3x3 matrix */
-void sa_Transpose3x3(double out[3][3], double in[3][3]) {
- int i, j;
- if (out != in) {
- for (i = 0; i < 3; i++)
- for (j = 0; j < 3; j++)
- out[i][j] = in[j][i];
- } else {
- double tt[3][3];
- for (i = 0; i < 3; i++)
- for (j = 0; j < 3; j++)
- tt[i][j] = in[j][i];
- for (i = 0; i < 3; i++)
- for (j = 0; j < 3; j++)
- out[i][j] = tt[i][j];
- }
-}
-
-/* Scale a 3 vector by the given ratio */
-void sa_Scale3(double out[3], double in[3], double rat) {
- out[0] = in[0] * rat;
- out[1] = in[1] * rat;
- out[2] = in[2] * rat;
-}
-
-/* Clamp a 3 vector to be +ve */
-void sa_Clamp3(double out[3], double in[3]) {
- int i;
- for (i = 0; i < 3; i++)
- out[i] = in[i] < 0.0 ? 0.0 : in[i];
-}
-
-/* Return the normal Delta E given two Lab values */
-double sa_LabDE(double *Lab0, double *Lab1) {
- double rv = 0.0, tt;
-
- tt = Lab0[0] - Lab1[0];
- rv += tt * tt;
- tt = Lab0[1] - Lab1[1];
- rv += tt * tt;
- tt = Lab0[2] - Lab1[2];
- rv += tt * tt;
-
- return sqrt(rv);
-}
-
-/* CIE XYZ to perceptual CIE 1976 L*a*b* */
-void
-sa_XYZ2Lab(icmXYZNumber *w, double *out, double *in) {
- double X = in[0], Y = in[1], Z = in[2];
- double x,y,z,fx,fy,fz;
-
- x = X/w->X;
- y = Y/w->Y;
- z = Z/w->Z;
-
- if (x > 0.008856451586)
- fx = pow(x,1.0/3.0);
- else
- fx = 7.787036979 * x + 16.0/116.0;
-
- if (y > 0.008856451586)
- fy = pow(y,1.0/3.0);
- else
- fy = 7.787036979 * y + 16.0/116.0;
-
- if (z > 0.008856451586)
- fz = pow(z,1.0/3.0);
- else
- fz = 7.787036979 * z + 16.0/116.0;
-
- out[0] = 116.0 * fy - 16.0;
- out[1] = 500.0 * (fx - fy);
- out[2] = 200.0 * (fy - fz);
-}
-
-/* - - - - - - - - - - - - - - - - - - - - - - - - */
-/* A sub-set of ludecomp code from numlib */
-
-int sa_lu_decomp(double **a, int n, int *pivx, double *rip) {
- int i, j;
- double *rscale, RSCALE[10];
-
- if (n <= 10)
- rscale = RSCALE;
- else
- rscale = dvector(0, n-1);
-
- for (i = 0; i < n; i++) {
- double big;
- for (big = 0.0, j=0; j < n; j++) {
- double temp;
- temp = fabs(a[i][j]);
- if (temp > big)
- big = temp;
- }
- if (fabs(big) <= DBL_MIN) {
- if (rscale != RSCALE)
- free_dvector(rscale, 0, n-1);
- return 1;
- }
- rscale[i] = 1.0/big;
- }
-
- for (*rip = 1.0, j = 0; j < n; j++) {
- double big;
- int k, bigi = 0;
-
- for (i = 0; i < j; i++) {
- double sum;
- sum = a[i][j];
- for (k = 0; k < i; k++)
- sum -= a[i][k] * a[k][j];
- a[i][j] = sum;
- }
-
- for (big = 0.0, i = j; i < n; i++) {
- double sum, temp;
- sum = a[i][j];
- for (k = 0; k < j; k++)
- sum -= a[i][k] * a[k][j];
- a[i][j] = sum;
- temp = rscale[i] * fabs(sum);
- if (temp >= big) {
- big = temp;
- bigi = i;
- }
- }
-
- if (j != bigi) {
- {
- double *temp;
- temp = a[bigi];
- a[bigi] = a[j];
- a[j] = temp;
- }
- *rip = -(*rip);
- rscale[bigi] = rscale[j];
- }
-
- pivx[j] = bigi;
- if (fabs(a[j][j]) <= DBL_MIN) {
- if (rscale != RSCALE)
- free_dvector(rscale, 0, n-1);
- return 1;
- }
-
- if (j != (n-1)) {
- double temp;
- temp = 1.0/a[j][j];
- for (i = j+1; i < n; i++)
- a[i][j] *= temp;
- }
- }
- if (rscale != RSCALE)
- free_dvector(rscale, 0, n-1);
- return 0;
-}
-
-void sa_lu_backsub(double **a, int n, int *pivx, double *b) {
- int i, j;
- int nvi;
-
- for (nvi = -1, i = 0; i < n; i++) {
- int px;
- double sum;
-
- px = pivx[i];
- sum = b[px];
- b[px] = b[i];
- if (nvi >= 0) {
- for (j = nvi; j < i; j++)
- sum -= a[i][j] * b[j];
- } else {
- if (sum != 0.0)
- nvi = i;
- }
- b[i] = sum;
- }
-
- for (i = (n-1); i >= 0; i--) {
- double sum;
- sum = b[i];
- for (j = i+1; j < n; j++)
- sum -= a[i][j] * b[j];
- b[i] = sum/a[i][i];
- }
-}
-
-int sa_lu_invert(double **a, int n) {
- int i, j;
- double rip;
- int *pivx, PIVX[10];
- double **y;
-
- if (n <= 10)
- pivx = PIVX;
- else
- pivx = ivector(0, n-1);
-
- if (sa_lu_decomp(a, n, pivx, &rip)) {
- if (pivx != PIVX)
- free_ivector(pivx, 0, n-1);
- return 1;
- }
-
- y = dmatrix(0, n-1, 0, n-1);
- for (i = 0; i < n; i++) {
- for (j = 0; j < n; j++) {
- y[i][j] = a[i][j];
- }
- }
-
- for (i = 0; i < n; i++) {
- for (j = 0; j < n; j++)
- a[i][j] = 0.0;
- a[i][i] = 1.0;
- sa_lu_backsub(y, n, pivx, a[i]);
- }
-
- free_dmatrix(y, 0, n-1, 0, n-1);
- if (pivx != PIVX)
- free_ivector(pivx, 0, n-1);
-
- return 0;
-}
-
-int sa_lu_psinvert(double **out, double **in, int m, int n) {
- int rv = 0;
- double **tr;
- double **sq;
-
- tr = dmatrix(0, n-1, 0, m-1);
- matrix_trans(tr, in, m, n);
-
- if (m > n) {
- sq = dmatrix(0, n-1, 0, n-1);
- if ((rv = matrix_mult(sq, n, n, tr, n, m, in, m, n)) == 0) {
- if ((rv = sa_lu_invert(sq, n)) == 0) {
- rv = matrix_mult(out, n, m, sq, n, n, tr, n, m);
- }
- }
- free_dmatrix(sq, 0, n-1, 0, n-1);
- } else {
- sq = dmatrix(0, m-1, 0, m-1);
- if ((rv = matrix_mult(sq, m, m, in, m, n, tr, n, m)) == 0) {
- if ((rv = sa_lu_invert(sq, m)) == 0) {
- rv = matrix_mult(out, n, m, tr, n, m, sq, m, m);
- }
- }
- free_dmatrix(sq, 0, m-1, 0, m-1);
- }
-
- free_dmatrix(tr, 0, n-1, 0, m-1);
- return rv;
-}
-
-
-#endif /* SALONEINSTLIB */
-/* ============================================================= */
+#endif /* UNIX_APPLE || NT */