Mercurial > hg > octave-nkf

diff scripts/linear-algebra/linsolve.m @ 20240:91e1da1d1918

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
linsolve.m: Calculate solution in a BLAS-optimized way (bug #44722). * linsolve.m: When transpose option is given, forward the transpose flag to BLAS rather than calculating the transpose in Octave and then forwarding the left-division operator to BLAS. Fixes inconsistency in BIST tests.
author Rik <rik@octave.org>
date Mon, 06 Apr 2015 04:18:47 -0700
parents 4197fc428c7d
children 014e942ac29f
line wrap: on
line diff
--- a/scripts/linear-algebra/linsolve.m
+++ b/scripts/linear-algebra/linsolve.m
@@ -87,28 +87,24 @@
     trans_A = false;
     if (isfield (opts, "TRANSA") && opts.TRANSA)
       trans_A = true;
-      A = A';
     endif
     if (isfield (opts, "POSDEF") && opts.POSDEF)
       A = matrix_type (A, "positive definite");
     endif
     if (isfield (opts, "LT") && opts.LT)
-      if (trans_A)
-        A = matrix_type (A, "upper");
-      else
-        A = matrix_type (A, "lower");
-      endif
-    endif
-    if (isfield (opts, "UT") && opts.UT)
-      if (trans_A)
-        A = matrix_type (A, "lower");
-      else
-        A = matrix_type (A, "upper");
-      endif
+      A = matrix_type (A, "lower");
+    elseif (isfield (opts, "UT") && opts.UT)
+      A = matrix_type (A, "upper");
     endif
   endif
 
-  x = A \ b;
+  ## This way is faster as the transpose is not calculated in Octave,
+  ## but forwarded as a flag option to BLAS.
+  if (trans_A)
+    x = A' \ b;
+  else
+    x = A \ b;
+  endif
 
   if (nargout > 1)
     if (issquare (A))
@@ -117,12 +113,13 @@
       R = 0;
     endif
   endif
+
 endfunction
 
 
 %!test
-%! n = 4;
-%! A = triu (rand (n));
+%! n = 10;
+%! A = triu (gallery ("condex", n));
 %! x = rand (n, 1);
 %! b = A' * x;
 %! opts.UT = true;