← Back to branch summary

~ubuntu-branches/ubuntu/edgy/ncbi-tools6/edgy

« back to all changes in this revision

Viewing changes to algo/blast/composition_adjustment/nlm_linear_algebra.c

  • Committer: Bazaar Package Importer
  • Author(s): Barry deFreese
  • Date: 2006-07-19 23:28:07 UTC
  • mfrom: (1.1.5 upstream)
  • Revision ID: james.westby@ubuntu.com-20060719232807-et3cdmcjgmnyleyx
Tags: 6.1.20060507-3ubuntu1
Re-merge with Debian

Show diffs side-by-side

added added

removed removed

Lines of Context:
algo/blast/composition_adjustment/nlm_linear_algebra.c
23
23
* ===========================================================================*/
24
24
 
25
25
/** @file nlm_linear_algebra.c
 
26
 * Basic matrix and vector operations
26
27
 *
27
28
 * @author E. Michael Gertz
28
 
 *
29
 
 * Basic matrix and vector operations
30
29
 */
31
30
#ifndef SKIP_DOXYGEN_PROCESSING
32
31
static char const rcsid[] =
33
 
    "$Id: nlm_linear_algebra.c,v 1.5 2005/12/01 13:49:43 gertz Exp $";
 
32
    "$Id: nlm_linear_algebra.c,v 1.7 2005/12/19 15:37:33 gertz Exp $";
34
33
#endif /* SKIP_DOXYGEN_PROCESSING */
35
34
 
36
35
#include <math.h>
38
37
#include <algo/blast/core/ncbi_std.h>
39
38
#include <algo/blast/composition_adjustment/nlm_linear_algebra.h>
40
39
 
41
 
/**
42
 
 * Create and return a new, dense matrix.  Elements of the matrix A
43
 
 * may be accessed as A[i][j]
44
 
 *
45
 
 * @param nrows     the number of rows for the new matrix.
46
 
 * @param ncols     the number of columns for the new matrix.
47
 
 */
 
40
 
 
41
/* Documented in nlm_linear_algebra.h. */
48
42
double **
49
43
Nlm_DenseMatrixNew(int nrows,
50
44
                   int ncols)
69
63
}
70
64
 
71
65
 
72
 
/**
73
 
 * Create and return a new, dense, lower-triangular matrix.  Elements
74
 
 * of the matrix A may be accessed as A[i][j] for i <= j.
75
 
 *
76
 
 * @param n         the dimension of the matrix.
77
 
 */
 
66
/* Documented in nlm_linear_algebra.h. */
78
67
double **
79
68
Nlm_LtriangMatrixNew(int n)
80
69
{
100
89
}
101
90
 
102
91
 
103
 
/**
104
 
 * Free a matrix created by Nlm_DenseMatrixNew or
105
 
 * Nlm_LtriangMatrixNew.
106
 
 *
107
 
 * @param mat       the matrix to be freed
108
 
 * @return          always NULL
109
 
 */
 
92
/* Documented in nlm_linear_algebra.h. */
110
93
void
111
94
Nlm_DenseMatrixFree(double *** mat)
112
95
{
118
101
}
119
102
 
120
103
 
121
 
/**
122
 
 * Create and return a new Int4 matrix.  Elements of the matrix A
123
 
 * may be accessed as A[i][j]
124
 
 *
125
 
 * @param nrows     the number of rows for the new matrix.
126
 
 * @param ncols     the number of columns for the new matrix.
127
 
 */
 
104
/* Documented in nlm_linear_algebra.h. */
128
105
Int4 ** Nlm_Int4MatrixNew(int nrows, int ncols)
129
106
{
130
107
    int i;             /* iteration index */
147
124
}
148
125
 
149
126
 
150
 
/**
151
 
 * Free a matrix created by Nlm_DenseMatrixNew or
152
 
 * Nlm_LtriangMatrixNew.
153
 
 *
154
 
 * @param mat       the matrix to be freed
155
 
 * @return          always NULL
156
 
 */
 
127
/* Documented in nlm_linear_algebra.h. */
157
128
void
158
129
Nlm_Int4MatrixFree(Int4 *** mat)
159
130
{
164
135
    *mat = NULL;
165
136
}
166
137
 
167
 
/**
168
 
 * Accessing only the lower triangular elements of the symmetric,
169
 
 * positive definite matrix A, compute a lower triangular matrix L
170
 
 * such that A = L L^T (Cholesky factorization.)  Overwrite the lower
171
 
 * triangle of A with L.
172
 
 *
173
 
 * This routine may be used with the Nlm_SolveLtriangPosDef routine to
174
 
 * solve systems of equations.
175
 
 *
176
 
 * @param A         the lower triangle of a symmetric, positive-definite
177
 
 *                  matrix
178
 
 * @param n         the size of A
179
 
 */
 
138
 
 
139
/* Documented in nlm_linear_algebra.h. */
180
140
void
181
141
Nlm_FactorLtriangPosDef(double ** A, int n)
182
142
{
201
161
}
202
162
 
203
163
 
204
 
/**
205
 
 * Solve the linear system L L\T y = b, where L is a non-singular
206
 
 * lower triangular matrix, usually computed using
207
 
 * the Nlm_FactorLtriangPosDef routine.
208
 
 *
209
 
 * @param x         on entry, the right hand size of the linear system
210
 
 *                  L L\T y = b; on exit the solution
211
 
 * @param n         the size of x
212
 
 * @param L         a non-singular lower triangular matrix
213
 
 */
214
 
void Nlm_SolveLtriangPosDef(double * x, int n,
 
164
/* Documented in nlm_linear_algebra.h. */
 
165
void Nlm_SolveLtriangPosDef(double x[], int n,
215
166
                            double ** L )
216
167
{
217
168
    int i, j;                   /* iteration indices */
239
190
}
240
191
 
241
192
 
242
 
/**
243
 
 * Compute the Euclidean norm (2-norm) of a vector.
244
 
 *
245
 
 * This routine is based on the (freely available) BLAS routine dnrm2,
246
 
 * which handles the scale of the elements of v in a stable fashion.
247
 
 *
248
 
 * @param v      a vector
249
 
 * @param n      the length of v
250
 
 */
 
193
/* Documented in nlm_linear_algebra.h. */
251
194
double
252
 
Nlm_EuclideanNorm(const double * v, int n)
 
195
Nlm_EuclideanNorm(const double v[], int n)
253
196
{
254
197
    double sum   = 1.0;   /* sum of squares of elements in v */
255
198
    double scale = 0.0;   /* a scale factor for the elements in v */
270
213
}
271
214
 
272
215
 
273
 
/**
274
 
 * Let y = y + alpha * x
275
 
 *
276
 
 * @param y         a vector
277
 
 * @param x         another vector
278
 
 * @param n         the length of x and y
279
 
 */
280
 
void Nlm_AddVectors(double * y, int n, double alpha, const double * x)
 
216
/* Documented in nlm_linear_algebra.h. */
 
217
void Nlm_AddVectors(double y[], int n, double alpha, const double x[])
281
218
{
282
219
    int i;                     /* iteration index */
283
220
 
287
224
}
288
225
 
289
226
 
290
 
/**
291
 
 * Given a nonnegative vector x and a nonnegative scalar max, returns
292
 
 * the largest value in [0, max] for which x + alpha * step_x >= 0.
293
 
 *
294
 
 * @param x         a vector with nonnegative elements
295
 
 * @param step_x    another vector
296
 
 * @param n         the size of x and step_x
297
 
 * @param max       a nonnegative scalar
298
 
 */
 
227
/* Documented in nlm_linear_algebra.h. */
299
228
double
300
 
Nlm_StepBound(const double * x, int n, const double * step_x, double max)
 
229
Nlm_StepBound(const double x[], int n, const double step_x[], double max)
301
230
{
302
231
    int i;                 /* iteration index */
303
232
    double alpha = max;    /* current largest permitted step */