1 /* ===========================================================================
2 *
3 * PUBLIC DOMAIN NOTICE
4 * National Center for Biotechnology Information
5 *
6 * This software/database is a "United States Government Work" under the
7 * terms of the United States Copyright Act. It was written as part of
8 * the author's official duties as a United States Government employee and
9 * thus cannot be copyrighted. This software/database is freely available
10 * to the public for use. The National Library of Medicine and the U.S.
11 * Government have not placed any restriction on its use or reproduction.
12 *
13 * Although all reasonable efforts have been taken to ensure the accuracy
14 * and reliability of the software and data, the NLM and the U.S.
15 * Government do not and cannot warrant the performance or results that
16 * may be obtained by using this software or data. The NLM and the U.S.
17 * Government disclaim all warranties, express or implied, including
18 * warranties of performance, merchantability or fitness for any particular
19 * purpose.
20 *
21 * Please cite the author in any work or product based on this material.
22 *
23 * ===========================================================================*/
24
25 /** @file nlm_linear_algebra.c
26 * Basic matrix and vector operations
27 *
28 * @author E. Michael Gertz
29 */
30
31 #include <math.h>
32 #include <stdlib.h>
33 #include <algo/blast/core/ncbi_std.h>
34 #include <algo/blast/composition_adjustment/nlm_linear_algebra.h>
35
36
37 /* Documented in nlm_linear_algebra.h. */
38 double **
Nlm_DenseMatrixNew(int nrows,int ncols)39 Nlm_DenseMatrixNew(int nrows,
40 int ncols)
41 {
42 int i; /* iteration index */
43 double ** mat; /* the new matrix */
44
45 mat = (double **) calloc(nrows, sizeof(double *));
46 if (mat != NULL) {
47 mat[0] = (double *) malloc((size_t) nrows *
48 (size_t) ncols * sizeof(double));
49 if (mat[0] != NULL) {
50 for (i = 1; i < nrows; i++) {
51 mat[i] = &mat[0][i * ncols];
52 }
53 } else {
54 free(mat);
55 mat = NULL;
56 }
57 }
58 return mat;
59 }
60
61
62 /* Documented in nlm_linear_algebra.h. */
63 double **
Nlm_LtriangMatrixNew(int n)64 Nlm_LtriangMatrixNew(int n)
65 {
66 int i; /* iteration index */
67 double ** L; /* the new, lower triangular matrix */
68 size_t nelts; /* the number of elements in
69 the matrix */
70 nelts = ((size_t) n * (n + 1))/2;
71
72 L = (double**) calloc(n, sizeof(double *));
73 if (L != NULL) {
74 L[0] = (double*) malloc(nelts * sizeof(double));
75 if (L[0] != NULL) {
76 for (i = 1; i < n; i++) {
77 L[i] = L[i - 1] + i;
78 }
79 } else {
80 free(L);
81 L = NULL;
82 }
83 }
84 return L;
85 }
86
87
88 /* Documented in nlm_linear_algebra.h. */
89 void
Nlm_DenseMatrixFree(double *** mat)90 Nlm_DenseMatrixFree(double *** mat)
91 {
92 if(*mat != NULL) {
93 free((*mat)[0]);
94 free(*mat);
95 }
96 *mat = NULL;
97 }
98
99
100 /* Documented in nlm_linear_algebra.h. */
Nlm_Int4MatrixNew(int nrows,int ncols)101 int ** Nlm_Int4MatrixNew(int nrows, int ncols)
102 {
103 int i; /* iteration index */
104 int ** mat; /* the new matrix */
105
106 mat = (int **) calloc(nrows, sizeof(int *));
107 if (mat != NULL) {
108 mat[0] = (int *) malloc((size_t) nrows *
109 (size_t) ncols * sizeof(int));
110 if (mat[0] != NULL) {
111 for (i = 1; i < nrows; i++) {
112 mat[i] = &mat[0][i * ncols];
113 }
114 } else {
115 free(mat);
116 mat = NULL;
117 }
118 }
119 return mat;
120 }
121
122
123 /* Documented in nlm_linear_algebra.h. */
124 void
Nlm_Int4MatrixFree(int *** mat)125 Nlm_Int4MatrixFree(int *** mat)
126 {
127 if(*mat != NULL) {
128 free((*mat)[0]);
129 free(*mat);
130 }
131 *mat = NULL;
132 }
133
134
135 /* Documented in nlm_linear_algebra.h. */
136 void
Nlm_FactorLtriangPosDef(double ** A,int n)137 Nlm_FactorLtriangPosDef(double ** A, int n)
138 {
139 int i, j, k; /* iteration indices */
140 double temp; /* temporary variable for intermediate
141 values in a computation */
142
143 for (i = 0; i < n; i++) {
144 for (j = 0; j < i; j++) {
145 temp = A[i][j];
146 for (k = 0; k < j; k++) {
147 temp -= A[i][k] * A[j][k];
148 }
149 A[i][j] = temp/A[j][j];
150 }
151 temp = A[i][i];
152 for (k = 0; k < i; k++) {
153 temp -= A[i][k] * A[i][k];
154 }
155 A[i][i] = sqrt(temp);
156 }
157 }
158
159
160 /* Documented in nlm_linear_algebra.h. */
Nlm_SolveLtriangPosDef(double x[],int n,double ** L)161 void Nlm_SolveLtriangPosDef(double x[], int n,
162 double ** L )
163 {
164 int i, j; /* iteration indices */
165 double temp; /* temporary variable for intermediate
166 values in a computation */
167
168 /* At point x = b in the equation L L\T y = b */
169
170 /* Forward solve; L z = b */
171 for (i = 0; i < n; i++) {
172 temp = x[i];
173 for (j = 0; j < i; j++) {
174 temp -= L[i][j] * x[j];
175 }
176 x[i] = temp/L[i][i];
177 }
178 /* Now x = z. Back solve the system L\T y = z */
179 for (j = n - 1; j >= 0; j--) {
180 x[j] /= L[j][j];
181 for (i = 0; i < j; i++) {
182 x[i] -= L[j][i] * x[j];
183 }
184 }
185 /* Now x = y, the solution to L L\T y = b */
186 }
187
188
189 /* Documented in nlm_linear_algebra.h. */
190 double
Nlm_EuclideanNorm(const double v[],int n)191 Nlm_EuclideanNorm(const double v[], int n)
192 {
193 double sum = 1.0; /* sum of squares of elements in v */
194 double scale = 0.0; /* a scale factor for the elements in v */
195 int i; /* iteration index */
196
197 for (i = 0; i < n; i++) {
198 if (v[i] != 0.0) {
199 double absvi = fabs(v[i]);
200 if (scale < absvi) {
201 sum = 1.0 + sum * (scale/absvi) * (scale/absvi);
202 scale = absvi;
203 } else {
204 sum += (absvi/scale) * (absvi/scale);
205 }
206 }
207 }
208 return scale * sqrt(sum);
209 }
210
211
212 /* Documented in nlm_linear_algebra.h. */
Nlm_AddVectors(double y[],int n,double alpha,const double x[])213 void Nlm_AddVectors(double y[], int n, double alpha, const double x[])
214 {
215 int i; /* iteration index */
216
217 for (i = 0; i < n; i++) {
218 y[i] += alpha * x[i];
219 }
220 }
221
222
223 /* Documented in nlm_linear_algebra.h. */
224 double
Nlm_StepBound(const double x[],int n,const double step_x[],double max)225 Nlm_StepBound(const double x[], int n, const double step_x[], double max)
226 {
227 int i; /* iteration index */
228 double alpha = max; /* current largest permitted step */
229
230 for (i = 0; i < n; i++) {
231 double alpha_i; /* a step to the boundary for the current i */
232
233 alpha_i = -x[i] / step_x[i];
234 if (alpha_i >= 0 && alpha_i < alpha) {
235 alpha = alpha_i;
236 }
237 }
238 return alpha;
239 }
240