00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include <string.h>
00038 #include <stdlib.h>
00039
00040 #ifdef HAVE_CONFIG_H
00041 #include "config.h"
00042 #endif
00043
00044 #include "clapack_lite.h"
00045 #include "matrix.h"
00046 #include "err.h"
00047 #include "ckd_alloc.h"
00048
00049 #ifndef WITH_LAPACK
00050 float64
00051 determinant(float32 **a, int32 n)
00052 {
00053 E_FATAL("No LAPACK library available, cannot compute determinant (FIXME)\n");
00054 return 0.0;
00055 }
00056 int32
00057 invert(float32 **ainv, float32 **a, int32 n)
00058 {
00059 E_FATAL("No LAPACK library available, cannot compute matrix inverse (FIXME)\n");
00060 return 0;
00061 }
00062 int32
00063 solve(float32 **a, float32 *b, float32 *out_x, int32 n)
00064 {
00065 E_FATAL("No LAPACK library available, cannot solve linear equations (FIXME)\n");
00066 return 0;
00067 }
00068
00069 void
00070 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
00071 {
00072 int32 i, j, k;
00073
00074 memset(c[0], 0, n*n*sizeof(float32));
00075 for (i = 0; i < n; ++i) {
00076 for (j = 0; j < n; ++j) {
00077 for (k = 0; k < n; ++k) {
00078 c[i][k] += a[i][j] * b[j][k];
00079 }
00080 }
00081 }
00082 }
00083 #else
00084
00085 float64
00086 determinant(float32 ** a, int32 n)
00087 {
00088 float32 **tmp_a;
00089 float64 det;
00090 char uplo;
00091 int32 info, i;
00092
00093
00094
00095
00096 tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
00097 memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
00098
00099 uplo = 'L';
00100 spotrf_(&uplo, &n, tmp_a[0], &n, &info);
00101 det = tmp_a[0][0];
00102
00103 for (i = 1; i < n; ++i)
00104 det *= tmp_a[i][i];
00105 ckd_free_2d((void **)tmp_a);
00106 if (info > 0)
00107 return -1.0;
00108 else
00109 return det * det;
00110 }
00111
00112 int32
00113 solve(float32 **a,
00114 float32 *b,
00115 float32 *out_x,
00116 int32 n)
00117 {
00118 char uplo;
00119 float32 **tmp_a;
00120 int32 info, nrhs;
00121
00122
00123
00124
00125 tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
00126 memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
00127 memcpy(out_x, b, n*sizeof(float32));
00128 uplo = 'L';
00129 nrhs = 1;
00130 sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, out_x, &n, &info);
00131 ckd_free_2d((void **)tmp_a);
00132
00133 if (info != 0)
00134 return -1;
00135 else
00136 return info;
00137 }
00138
00139
00140 int32
00141 invert(float32 ** ainv, float32 ** a, int32 n)
00142 {
00143 char uplo;
00144 float32 **tmp_a;
00145 int32 info, nrhs, i;
00146
00147
00148 memset(ainv[0], 0, sizeof(float32) * n * n);
00149 for (i = 0; i < n; i++)
00150 ainv[i][i] = 1.0;
00151
00152
00153
00154 tmp_a = (float32 **)ckd_calloc_2d(n, n, sizeof(float32));
00155 memcpy(tmp_a[0], a[0], n*n*sizeof(float32));
00156 uplo = 'L';
00157 nrhs = n;
00158 sposv_(&uplo, &n, &nrhs, tmp_a[0], &n, ainv[0], &n, &info);
00159 ckd_free_2d((void **)tmp_a);
00160
00161 if (info != 0)
00162 return -1;
00163 else
00164 return info;
00165 }
00166
00167 void
00168 matrixmultiply(float32 ** c, float32 ** a, float32 ** b, int32 n)
00169 {
00170 char side, uplo;
00171 float32 alpha;
00172
00173 side = 'L';
00174 uplo = 'L';
00175 alpha = 1.0;
00176 ssymm_(&side, &uplo, &n, &n, &alpha, a[0], &n, b[0], &n, &alpha, c[0], &n);
00177 }
00178
00179 #endif
00180
00181 void
00182 outerproduct(float32 ** a, float32 * x, float32 * y, int32 len)
00183 {
00184 int32 i, j;
00185
00186 for (i = 0; i < len; ++i) {
00187 a[i][i] = x[i] * y[i];
00188 for (j = i + 1; j < len; ++j) {
00189 a[i][j] = x[i] * y[j];
00190 a[j][i] = x[j] * y[i];
00191 }
00192 }
00193 }
00194
00195 void
00196 scalarmultiply(float32 ** a, float32 x, int32 n)
00197 {
00198 int32 i, j;
00199
00200 for (i = 0; i < n; ++i) {
00201 a[i][i] *= x;
00202 for (j = i+1; j < n; ++j) {
00203 a[i][j] *= x;
00204 a[j][i] *= x;
00205 }
00206 }
00207 }
00208
00209 void
00210 matrixadd(float32 ** a, float32 ** b, int32 n)
00211 {
00212 int32 i, j;
00213
00214 for (i = 0; i < n; ++i)
00215 for (j = 0; j < n; ++j)
00216 a[i][j] += b[i][j];
00217 }
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240