2 * Introduction to lapack::
3 * Functions and Variables for lapack::
6 @node Introduction to lapack, Functions and Variables for lapack, Package lapack, Package lapack
7 @section Introduction to lapack
9 @code{lapack} is a Common Lisp translation (via the program @code{f2cl}) of the Fortran library LAPACK,
10 as obtained from the SLATEC project.
12 @opencatbox{Categories:}
13 @category{Numerical methods}
14 @category{Share packages}
15 @category{Package lapack}
19 @node Functions and Variables for lapack, , Introduction to lapack, Package lapack
20 @section Functions and Variables for lapack
23 @deffn {Function} dgeev @
24 @fname{dgeev} (@var{A}) @
25 @fname{dgeev} (@var{A}, @var{right_p}, @var{left_p})
27 Computes the eigenvalues and, optionally, the eigenvectors of a matrix @var{A}.
28 All elements of @var{A} must be integer or floating point numbers.
29 @var{A} must be square (same number of rows and columns).
30 @var{A} might or might not be symmetric.
32 @code{dgeev(@var{A})} computes only the eigenvalues of @var{A}.
33 @code{dgeev(@var{A}, @var{right_p}, @var{left_p})} computes the eigenvalues of @var{A}
34 and the right eigenvectors when @math{@var{right_p} = @code{true}}
35 and the left eigenvectors when @math{@var{left_p} = @code{true}}.
37 A list of three items is returned.
38 The first item is a list of the eigenvalues.
39 The second item is @code{false} or the matrix of right eigenvectors.
40 The third item is @code{false} or the matrix of left eigenvectors.
42 The right eigenvector @math{v(j)} (the @math{j}-th column of the right eigenvector matrix) satisfies
44 @math{A . v(j) = lambda(j) . v(j)}
46 where @math{lambda(j)} is the corresponding eigenvalue.
47 The left eigenvector @math{u(j)} (the @math{j}-th column of the left eigenvector matrix) satisfies
49 @math{u(j)**H . A = lambda(j) . u(j)**H}
51 where @math{u(j)**H} denotes the conjugate transpose of @math{u(j)}.
52 The Maxima function @code{ctranspose} computes the conjugate transpose.
54 The computed eigenvectors are normalized to have Euclidean norm
55 equal to 1, and largest component has imaginary part equal to zero.
62 @c M : matrix ([9.5, 1.75], [3.25, 10.45]);
64 @c [L, v, u] : dgeev (M, true, true);
65 @c D : apply (diag_matrix, L);
67 @c transpose (u) . M - D . transpose (u);
70 (%i1) load ("lapack")$
71 (%i2) fpprintprec : 6;
73 (%i3) M : matrix ([9.5, 1.75], [3.25, 10.45]);
78 (%o4) [[7.54331, 12.4067], false, false]
79 (%i5) [L, v, u] : dgeev (M, true, true);
80 [ - .666642 - .515792 ]
81 (%o5) [[7.54331, 12.4067], [ ],
83 [ - .856714 - .745378 ]
86 (%i6) D : apply (diag_matrix, L);
94 (%i8) transpose (u) . M - D . transpose (u);
100 @opencatbox{Categories:}
101 @category{Package lapack}
107 @deffn {Function} dgeqrf (@var{A})
109 Computes the QR decomposition of the matrix @var{A}.
110 All elements of @var{A} must be integer or floating point numbers.
111 @var{A} may or may not have the same number of rows and columns.
113 A list of two items is returned.
114 The first item is the matrix @var{Q}, which is a square, orthonormal matrix
115 which has the same number of rows as @var{A}.
116 The second item is the matrix @var{R}, which is the same size as @var{A},
117 and which has all elements equal to zero below the diagonal.
118 The product @code{@var{Q} . @var{R}}, where "." is the noncommutative multiplication operator,
119 is equal to @var{A} (ignoring floating point round-off errors).
124 @c M : matrix ([1, -3.2, 8], [-11, 2.7, 5.9]);
125 @c [q, r] : dgeqrf (M);
130 (%i1) load ("lapack") $
131 (%i2) fpprintprec : 6 $
132 (%i3) M : matrix ([1, -3.2, 8], [-11, 2.7, 5.9]) $
133 (%i4) [q, r] : dgeqrf (M);
134 [ - .0905357 .995893 ]
137 [ - 11.0454 2.97863 5.15148 ]
139 [ 0 - 2.94241 8.50131 ]
141 [ - 7.77156E-16 1.77636E-15 - 8.88178E-16 ]
143 [ 0.0 - 1.33227E-15 8.88178E-16 ]
144 (%i6) mat_norm (%, 1);
148 @opencatbox{Categories:}
149 @category{Package lapack}
155 @deffn {Function} dgesv (@var{A}, @var{b})
157 Computes the solution @var{x} of the linear equation @math{@var{A} @var{x} = @var{b}},
158 where @var{A} is a square matrix, and @var{b} is a matrix of the same number of rows
159 as @var{A} and any number of columns.
160 The return value @var{x} is the same size as @var{b}.
162 The elements of @var{A} and @var{b} must evaluate to real floating point numbers via @code{float};
163 thus elements may be any numeric type, symbolic numerical constants, or expressions which evaluate to floats.
164 The elements of @var{x} are always floating point numbers.
165 All arithmetic is carried out as floating point operations.
167 @code{dgesv} computes the solution via the LU decomposition of @var{A}.
171 @code{dgesv} computes the solution of the linear equation @math{@var{A} @var{x} = @var{b}}.
174 @c A : matrix ([1, -2.5], [0.375, 5]);
175 @c b : matrix ([1.75], [-0.625]);
177 @c dlange (inf_norm, b - A . x);
180 (%i1) A : matrix ([1, -2.5], [0.375, 5]);
184 (%i2) b : matrix ([1.75], [-0.625]);
188 (%i3) x : dgesv (A, b);
189 [ 1.210526315789474 ]
191 [ - 0.215789473684211 ]
192 (%i4) dlange (inf_norm, b - A.x);
196 @var{b} is a matrix with the same number of rows as @var{A} and any number of columns.
197 @var{x} is the same size as @var{b}.
200 @c A : matrix ([1, -0.15], [1.82, 2]);
201 @c b : matrix ([3.7, 1, 8], [-2.3, 5, -3.9]);
203 @c dlange (inf_norm, b - A . x);
206 (%i1) A : matrix ([1, -0.15], [1.82, 2]);
210 (%i2) b : matrix ([3.7, 1, 8], [-2.3, 5, -3.9]);
214 (%i3) x : dgesv (A, b);
215 [ 3.103827540695117 1.20985481742191 6.781786185657722 ]
217 [ -3.974483062032557 1.399032116146062 -8.121425428948527 ]
218 (%i4) dlange (inf_norm, b - A . x);
219 (%o4) 1.1102230246251565E-15
222 The elements of @var{A} and @var{b} must evaluate to real floating point numbers.
225 @c A : matrix ([5, -%pi], [1b0, 11/17]);
226 @c b : matrix ([%e], [sin(1)]);
228 @c dlange (inf_norm, b - A . x);
231 (%i1) A : matrix ([5, -%pi], [1b0, 11/17]);
237 (%i2) b : matrix ([%e], [sin(1)]);
241 (%i3) x : dgesv (A, b);
242 [ 0.690375643155986 ]
244 [ 0.233510982552952 ]
245 (%i4) dlange (inf_norm, b - A . x);
246 (%o4) 2.220446049250313E-16
249 @opencatbox{Categories:}
250 @category{Package lapack}
251 @category{Linear equations}
257 @deffn {Function} dgesvd @
258 @fname{dgesvd} (@var{A}) @
259 @fname{dgesvd} (@var{A}, @var{left_p}, @var{right_p})
261 Computes the singular value decomposition (SVD) of a matrix @var{A},
262 comprising the singular values and, optionally, the left and right singular vectors.
263 All elements of @var{A} must be integer or floating point numbers.
264 @var{A} might or might not be square (same number of rows and columns).
266 Let @math{m} be the number of rows, and @math{n} the number of columns of @var{A}.
267 The singular value decomposition of @var{A} comprises three matrices,
268 @var{U}, @var{Sigma}, and @var{V^T},
271 @c this code breaks texi2pdf: @math{@var{A} = @var{U} . @var{Sigma} . @var{V^T}}
272 @math{@var{A} = @var{U} . @var{Sigma} . @var{V}^T}
274 where @var{U} is an @math{m}-by-@math{m} unitary matrix,
275 @var{Sigma} is an @math{m}-by-@math{n} diagonal matrix,
276 and @var{V^T} is an @math{n}-by-@math{n} unitary matrix.
278 Let @math{sigma[i]} be a diagonal element of @math{Sigma},
279 that is, @math{@var{Sigma}[i, i] = @var{sigma}[i]}.
280 The elements @math{sigma[i]} are the so-called singular values of @var{A};
281 these are real and nonnegative, and returned in descending order.
282 The first @math{min(m, n)} columns of @var{U} and @var{V} are
283 the left and right singular vectors of @var{A}.
284 Note that @code{dgesvd} returns the transpose of @var{V}, not @var{V} itself.
286 @code{dgesvd(@var{A})} computes only the singular values of @var{A}.
287 @code{dgesvd(@var{A}, @var{left_p}, @var{right_p})} computes the singular values of @var{A}
288 and the left singular vectors when @math{@var{left_p} = @code{true}}
289 and the right singular vectors when @math{@var{right_p} = @code{true}}.
291 A list of three items is returned.
292 The first item is a list of the singular values.
293 The second item is @code{false} or the matrix of left singular vectors.
294 The third item is @code{false} or the matrix of right singular vectors.
301 @c M: matrix([1, 2, 3], [3.5, 0.5, 8], [-1, 2, -3], [4, 9, 7]);
303 @c [sigma, U, VT] : dgesvd (M, true, true);
307 @c genmatrix(lambda ([i, j], if i=j then sigma[i] else 0),
309 @c U . Sigma . VT - M;
310 @c transpose (U) . U;
311 @c VT . transpose (VT);
314 (%i1) load ("lapack")$
315 (%i2) fpprintprec : 6;
317 (%i3) M: matrix([1, 2, 3], [3.5, 0.5, 8], [-1, 2, -3], [4, 9, 7]);
326 (%o4) [[14.4744, 6.38637, .452547], false, false]
327 (%i5) [sigma, U, VT] : dgesvd (M, true, true);
328 (%o5) [[14.4744, 6.38637, .452547],
329 [ - .256731 .00816168 .959029 - .119523 ]
331 [ - .526456 .672116 - .206236 - .478091 ]
333 [ .107997 - .532278 - .0708315 - 0.83666 ]
335 [ - .803287 - .514659 - .180867 .239046 ]
336 [ - .374486 - .538209 - .755044 ]
338 [ .130623 - .836799 0.5317 ]]
340 [ - .917986 .100488 .383672 ]
341 (%i6) m : length (U);
343 (%i7) n : length (VT);
346 genmatrix(lambda ([i, j], if i=j then sigma[i] else 0),
355 (%i9) U . Sigma . VT - M;
356 [ 1.11022E-15 0.0 1.77636E-15 ]
358 [ 1.33227E-15 1.66533E-15 0.0 ]
360 [ - 4.44089E-16 - 8.88178E-16 4.44089E-16 ]
362 [ 8.88178E-16 1.77636E-15 8.88178E-16 ]
363 (%i10) transpose (U) . U;
364 [ 1.0 5.55112E-17 2.498E-16 2.77556E-17 ]
366 [ 5.55112E-17 1.0 5.55112E-17 4.16334E-17 ]
368 [ 2.498E-16 5.55112E-17 1.0 - 2.08167E-16 ]
370 [ 2.77556E-17 4.16334E-17 - 2.08167E-16 1.0 ]
371 (%i11) VT . transpose (VT);
372 [ 1.0 0.0 - 5.55112E-17 ]
374 (%o11) [ 0.0 1.0 5.55112E-17 ]
376 [ - 5.55112E-17 5.55112E-17 1.0 ]
379 @opencatbox{Categories:}
380 @category{Package lapack}
387 @deffn {Function} dlange (@var{norm}, @var{A})
388 @deffnx {Function} zlange (@var{norm}, @var{A})
390 Computes a norm or norm-like function of the matrix @var{A}. If
391 @var{A} is a real matrix, use @code{dlange}. For a matrix with
392 complex elements, use @code{zlange}.
394 @code{norm} specifies the kind of norm to be computed:
397 Compute @math{max(abs(A(i, j)))} where @math{i} and @math{j} range over
398 the rows and columns, respectively, of @var{A}.
399 Note that this function is not a proper matrix norm.
401 Compute the @math{L[1]} norm of @var{A},
402 that is, the maximum of the sum of the absolute value of elements in each column.
404 Compute the @math{L[inf]} norm of @var{A},
405 that is, the maximum of the sum of the absolute value of elements in each row.
407 Compute the Frobenius norm of @var{A},
408 that is, the square root of the sum of squares of the matrix elements.
411 @opencatbox{Categories:}
412 @category{Package lapack}
418 @deffn {Function} dgemm @
419 @fname{dgemm} (@var{A}, @var{B}) @
420 @fname{dgemm} (@var{A}, @var{B}, @var{options})
421 Compute the product of two matrices and optionally add the product to
424 In the simplest form, @code{dgemm(@var{A}, @var{B})} computes the
425 product of the two real matrices, @var{A} and @var{B}.
427 In the second form, @code{dgemm} computes the @math{@var{alpha} *
428 @var{A} * @var{B} + @var{beta} * @var{C}} where @var{A}, @var{B},
429 @var{C} are real matrices of the appropriate sizes and @var{alpha} and
430 @var{beta} are real numbers. Optionally, @var{A} and/or @var{B} can
431 be transposed before computing the product. The extra parameters are
432 specified by optional keyword arguments: The keyword arguments are
433 optional and may be specified in any order. They all take the form
434 @code{key=val}. The keyword arguments are:
438 The matrix @var{C} that should be added. The default is @code{false},
439 which means no matrix is added.
441 The product of @var{A} and @var{B} is multiplied by this value. The
444 If a matrix @var{C} is given, this value multiplies @var{C} before it
445 is added. The default value is 0, which implies that @var{C} is not
446 added, even if @var{C} is given. Hence, be sure to specify a non-zero
447 value for @var{beta}.
449 If @code{true}, the transpose of @var{A} is used instead of @var{A}
450 for the product. The default is @code{false}.
452 If @code{true}, the transpose of @var{B} is used instead of @var{B}
453 for the product. The default is @code{false}.
458 @c A : matrix([1,2,3],[4,5,6],[7,8,9]);
459 @c B : matrix([-1,-2,-3],[-4,-5,-6],[-7,-8,-9]);
460 @c C : matrix([3,2,1],[6,5,4],[9,8,7]);
463 @c dgemm(A,B,transpose_a=true);
465 @c dgemm(A,B,c=C,beta=1);
467 @c dgemm(A,B,c=C,beta=1, alpha=-1);
471 (%i1) load ("lapack")$
472 (%i2) A : matrix([1,2,3],[4,5,6],[7,8,9]);
478 (%i3) B : matrix([-1,-2,-3],[-4,-5,-6],[-7,-8,-9]);
481 (%o3) [ - 4 - 5 - 6 ]
484 (%i4) C : matrix([3,2,1],[6,5,4],[9,8,7]);
491 [ - 30.0 - 36.0 - 42.0 ]
493 (%o5) [ - 66.0 - 81.0 - 96.0 ]
495 [ - 102.0 - 126.0 - 150.0 ]
499 (%o6) [ - 66 - 81 - 96 ]
501 [ - 102 - 126 - 150 ]
502 (%i7) dgemm(A,B,transpose_a=true);
503 [ - 66.0 - 78.0 - 90.0 ]
505 (%o7) [ - 78.0 - 93.0 - 108.0 ]
507 [ - 90.0 - 108.0 - 126.0 ]
508 (%i8) transpose(A) . B;
511 (%o8) [ - 78 - 93 - 108 ]
514 (%i9) dgemm(A,B,c=C,beta=1);
515 [ - 27.0 - 34.0 - 41.0 ]
517 (%o9) [ - 60.0 - 76.0 - 92.0 ]
519 [ - 93.0 - 118.0 - 143.0 ]
523 (%o10) [ - 60 - 76 - 92 ]
526 (%i11) dgemm(A,B,c=C,beta=1, alpha=-1);
529 (%o11) [ 72.0 86.0 100.0 ]
531 [ 111.0 134.0 157.0 ]
540 @opencatbox{Categories:}
541 @category{Package lapack}
547 @deffn {Function} zgeev @
548 @fname{zgeev} (@var{A}) @
549 @fname{zgeev} (@var{A}, @var{right_p}, @var{left_p})
551 Like @code{dgeev}, but the matrix @var{A} is complex.
553 @opencatbox{Categories:}
554 @category{Package lapack}
560 @deffn {Function} zheev @
561 @fname{zheev} (@var{A}) @
562 @fname{zheev} (@var{A}, @var{eigvec_p})
564 Like @code{dgeev}, but the matrix @var{A} is assumed to be a square
565 complex Hermitian matrix. If @var{eigvec_p} is @code{true}, then the
566 eigenvectors of the matrix are also computed.
568 No check is made that the matrix @var{A} is, in fact, Hermitian.
570 A list of two items is returned, as in @code{dgeev}: a list of
571 eigenvalues, and @code{false} or the matrix of the eigenvectors.
576 @c [9.14 +%i*0.00 , -4.37 -%i*9.22 , -1.98 -%i*1.72 , -8.96 -%i*9.50],
577 @c [-4.37 +%i*9.22 , -3.35 +%i*0.00 , 2.25 -%i*9.51 , 2.57 +%i*2.40],
578 @c [-1.98 +%i*1.72 , 2.25 +%i*9.51 , -4.82 +%i*0.00 , -3.24 +%i*2.04],
579 @c [-8.96 +%i*9.50 , 2.57 -%i*2.40 , -3.24 -%i*2.04 , 8.44 +%i*0.00]);
584 An example of computing the eigenvalues and then eigenvalues and
585 eigenvectors of an Hermitian matrix.
587 (%i1) load("lapack")$
589 [9.14 +%i*0.00 , -4.37 -%i*9.22 , -1.98 -%i*1.72 , -8.96 -%i*9.50],
590 [-4.37 +%i*9.22 , -3.35 +%i*0.00 , 2.25 -%i*9.51 , 2.57 +%i*2.40],
591 [-1.98 +%i*1.72 , 2.25 +%i*9.51 , -4.82 +%i*0.00 , -3.24 +%i*2.04],
592 [-8.96 +%i*9.50 , 2.57 -%i*2.40 , -3.24 -%i*2.04 , 8.44 +%i*0.00]);
594 [ 9.14 (- 9.22 %i) - 4.37 (- 1.72 %i) - 1.98 (- 9.5 %i) - 8.96 ]
596 [ 9.22 %i - 4.37 - 3.35 2.25 - 9.51 %i 2.4 %i + 2.57 ]
598 [ 1.72 %i - 1.98 9.51 %i + 2.25 - 4.82 2.04 %i - 3.24 ]
600 [ 9.5 %i - 8.96 2.57 - 2.4 %i (- 2.04 %i) - 3.24 8.44 ]
602 (%o3) [[- 16.00474647209473, - 6.764970154793324, 6.665711453507098,
603 25.51400517338097], false]
604 (%i4) E:zheev(A,true)$
606 (%o5) [- 16.00474647209474, - 6.764970154793325, 6.665711453507101,
609 [ 0.2674650533172745 %i + 0.2175453586665017 ]
611 [ 0.002696730886619885 %i + 0.6968836773391712 ]
613 [ (- 0.6082406376714117 %i) - 0.01210614292697931 ]
615 [ 0.1593081858095037 ]
616 [ 0.2644937470667444 %i + 0.4773693349937472 ]
618 [ (- 0.2852389036031621 %i) - 0.1414362742011673 ]
620 [ 0.2654607680986639 %i + 0.4467818117184174 ]
622 [ 0.5750762708542709 ]
623 [ 0.2810649767305922 %i - 0.1335263928245182 ]
625 [ 0.2866310132869556 %i - 0.4536971347853274 ]
627 [ (- 0.2933684323754295 %i) - 0.4954972425541057 ]
629 [ 0.5325337537576771 ]
630 [ (- 0.5737316575503476 %i) - 0.3966146799427706 ]
632 [ 0.01826502619021457 %i + 0.3530557704387017 ]
634 [ 0.1673700900085425 %i + 0.01476684746229564 ]
636 [ 0.6002632636961784 ]
639 @opencatbox{Categories:}
640 @category{Package lapack}