Writing a procedure for a matrix

[miktexuser]

How is this done?

If I use a procedure from the linalg package as a model fro writing a new procedure, I get an error when I try to use it.

For example, this is a renamed version of the linalg procedure.

symbolic procedure matrow(in_mat);
%
% Finds row dimension of a matrix.
%
begin
return first size_of_matrix(in_mat);
end;

If I try to use this in my program I get an error that matrow is not defined. I have loaded linalg so it should know what
first size_of_matrix is.

Does anyone know what the problem is?

Thanks.

[Rainer Schöpf]

This topic is covered in the "REDUCE Symbolic Mode Primer". An HTML version is linked on the REDUCE web site, and the TeX and PDF files can be found in the doc/primers subdirectory of the subversion repository.

Rainer

[Arthur Norman]

I will explain what is wrong for you, but the solution is going to be a
different matter!

Reduce has its language, but that language has two modes. One is referred
to as "algebraic" mode and is the one used by end-users where the
arguments you pass and the values you return from procedures are algebraic
formulae. The second is "symbolic" mode and is used wither by issuing a
directive "symbolic;" or by prefixing a statement or command with that
word. Symbolic mode is for Reduce developers and implementers and although
the syntax used is pretty much the same as algebraic mode the meaning is
quite different - the values you work with are the data-structures that
Reduce uses to represent formulae and the like. If you look in the
doc/primers part of the Reduce tree on Sourceforge you will find several
documents that aim to help those who wish to delve into Reduce internals
and chase bugs there, add new facilities or even just understand how
Reduce works.

So for MUCH of Reduce (eg much of the linalg package in particular) the
code you find is NOT user-level "algebraic mode" code that works at the
abstraction level of formulae, it is implementer-level "symbolic mode" and
works at the level of raw data structures.

There are a variety of schemes that arrange that capabilities implemented
in symbolic mode can be reflected and made avaialble to algebraic mode (ie
normal) users. After all that is the purpose of symbolic mode - to be the
implementation of all the algebra you are going to use.

In a hypothetical system other than Reduce the internal behaviour would
all have been coded up in some language that did not look at all like the
user-level one - for instance one could have an algebra system where the
implementation was in C or Java or rae Lisp or ... basically anything.
Reduce is perhaps unusual (and I view it as being in a good way!) that
(almost) the same syntax is used at both levels - but iof course it
represents an opportunity for confusion for you until you understand - and
others until they have read this posting!

Note that [trunk]/doc/manual/manual.pdf is the main Reduce manual. Of
course it is so huge that it is hard to take it all in at once!

In terms of what it looks as if you are trying to do in the particular
caswe you cite, I note
matrix a(3,7);
n := length a;
yields the list {3,7}, and then first n => 3 and second n => 7.

Arthur

On Wed, 21 Oct 2015, miktexuser wrote:

How is this done?

If I use a procedure from the linalg package as a model fro writing a new procedure, I get an error when I try to use it.

For example, this is a renamed version of the linalg procedure.

symbolic procedure matrow(in_mat);
%
% Finds row dimension of a matrix.
%
begin
return first size_of_matrix(in_mat);
end;

If I try to use this in my program I get an error that matrow is not defined. I have loaded linalg so it should know what
first size_of_matrix is.

Does anyone know what the problem is?

Thanks.

---

Writing a procedure for a matrix

---

Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/reduce-algebra/discussion/899364/

To unsubscribe from further messages, please visit https://sourceforge.net/auth/subscriptions/

[miktexuser]

Arthur

Thanks for the explanation. I was doing ok with normal procedures. The problem came when I wanted to use matrices in the produres which are not treated as in normal programs.

2 questions for symbolic procedures if I may.

If I have a matrix a, what's the best way to get the i,j term and set it to a variable eg. I want x = a(i,j).

If I have two matrices a and b in a procedure and use them to calculate another matrix c, how is c calculated and returned as a matrix?

Eg for all terms in a and b, find c(i,j) = a(i,j) * b(i,j), then return c or find a(i,j) = a(i,j) * b(i,j) and return the new a. Of course a and b are the same size in this example.

Thank you.

[Raffaele Vitolo]

On 22/10/15 04:24, miktexuser wrote:

Norman:

Thanks for the explanation. I was doing ok with normal procedures. The problem came when I wanted to use matrices in the produres which are not treated as in normal programs.

As I understand your question, a possible answer might be found in the
book on Reduce by MacCallum-Wright, they define the following procedures:
symbolic operator getmatel;
symbolic procedure getmatel(mat_d,i,j);
% MacCallum-Wright p.245
% handling of local matrices
nth(nth(cdr mat_d,i),j);

symbolic operator setmatel;
symbolic procedure setmatel(mat_d,i,j,val);
% MacCallum-Wright p. 245
% handling of local matrices
<<

rplaca(pnth(nth(cdr mat_d,i),j),val);
mat_d

;

The above procedures can be used inside other algebraic or symbolic
procedures to get and set elements of a matrix. Note the use of
`symbolic operator getmatel' that let you use a symbolic procedure in
algebraic mode too.

Raf

2 questions for symbolic procedures if I may.

If I have a matrix a, what's the best way to get the i,j term and set it to a variable eg. I want x = a(i,j).

If I have two matrices a and b in a procedure and use them to calculate another matrix c, how is c calculated and returned as a matrix?

Eg for all terms in a and b, find c(i,j) = a(i,j) * b(i,j), then return c or find a(i,j) = a(i,j) * b(i,j) and return the new a.

Thank you.

---

Writing a procedure for a matrix

---

Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/reduce-algebra/discussion/899364/

To unsubscribe from further messages, please visit https://sourceforge.net/auth/subscriptions/

[Arthur Norman]

On Thu, 22 Oct 2015, miktexuser wrote:

Arthur Norman:
Thanks for the explanation. I was doing ok with normal procedures. The
problem came when I wanted to use matrices in the produres which are not
treated as in normal programs.

2 questions for symbolic procedures if I may.

Raf has pointed you in the right direction for this, but use of SYMBOLIC
procedures should be viewed as something for fairly experienced users and
developers - eg those who have read all the primers and are happy coping
with the idea of the interface between user and system-builder mode
coding.

I had - perhaps incorrectly - sort of believed that yuou were just
starting out with Reduce in which case I would really suggest that you
leave symbolic mode stricyly alone until you HAVE to use it. And that your
first goes at solving your problems use algebraic mode Reduce somewhat as
a scripting language where you prepare a script to solve your particular
problem in a file and read that in... in such cases you often do not even
reach the state of needing to write "algebraic mode" procedures - you just
transcribe out what you want done in each case you want it done.

That strategy may let you get involved in getting the computations that
are relevant to you done with the least effort and pain - and then
wrapping the techniques up to produce a nice body of procedures and
interfaces and so on that would make what you have done easy for others to
take advantage of can come just a bit later.

I may have grossly misunderstood your level of experience with computer
algebra in general and Reduce in particular, in which case apologies!

          Arthur

[miktexuser]

Thank you Raf and Arthur.

Raf - The two routines are what I needed to access matrix elements in procedures.

Arthur - I have used procedures in Maple where it is easy to use matrices in procedures. I didn't realise that matrices are accessed differently in Reduce procedures.

Thanks to you both for your help.

[miktexuser]

How do I access pi = 3.142.. in a symbolic routine so that it knows the value of pi inside the routine (the same issue would arise for other reserved identifiers). I could pass its value as a parameter or calculate its value inside the routine, but is there a better way?

Thanks.

[Arthur Norman]

On Fri, 30 Oct 2015, miktexuser wrote:

How do I access pi = 3.142.. in a symbolic routine so that it knows the
value of pi inside the routine (the same issue would arise for other
reserved identifiers). I could pass its value as a parameter or calculate
its value inside the routine, but is there a better way?

Thanks.

You specifically say a SYMBOLIC procedure so you are working at the
implementation level in terms of raw data-structures. In that case perhaps
the easiest answer is that you can "simp 'pi" because in SYMBOLIC mode
simp will take any prefix form and evaluate it to get a standard quotient.
Obviously the value that you get will then depend on whether you have set
"in rounded" or what precision you are working to, so depending on those
sorts of things you would get back one of
((((pi . 1) . 1)) . 1)
((!:rd!: . 3.14159) . 1)
((!:rd!: 118686014801135580649625 . -75) . 1)
etc

packages/arith/bfelem.red is where the implementatiomn of lots of the
elementary constants can be found.

Is this what you wanted?
Arthur

[miktexuser]

Thank you Arthur. That is what I wanted.

[Thomas Sturm]

1: procedure hugo(x);
1: evalf(pi) + x;

hugo

2: hugo(0);

 27633741218861 
――――――――――――――――
 8796093022208  

3: hugo(1);

 36429834241069 
――――――――――――――――
 8796093022208  

4: precision 100;

12

5: hugo(0);

 4565260936322919454209049117203702099241912772600765 
――――――――――――――――――――――――――――――――――――――――――――――――――――――
 1453167689040254123888391366505616427928095106224401 

Best,

Thomas

[miktexuser]

Thank you Thomas. My question wasn't clear. If I use your method as follows.

symbolic operator hugo;
symbolic procedure hugo(x);
evalf(pi) + x;

hugo (2);

This gives an error that pi is undefined in hugo.

So my question was if there is a way of avoiding this without passing pi as a parameter.

Thanks.

[Thomas Sturm]

If you do not want to switch to rounded but use rational approximations of pi (subject to "precision" in your symbolic code, then you can use the symbolic code underlying evalf, like

23* simp evalf0 '(pi);

(27633741218861 . 8796093022208)

[miktexuser]

Thank you Thomas.

[miktexuser]

I have been reading McCallum/Wrights chapter 9 and trying the examples. The matrix routines work except norm1.

Here is the test code using their codes.

load_package matrix;

symbolic operator GetMatEl;
symbolic procedure GetMatEl(M, i, j);
% Returns M(i, j) where M is a LOCAL matrix.
% Based on getmatelem in matr.red
% Internally , M =(mat (m11 m12 ... )(m21 m22 ... ) ... )
<<
LocMatCheck(M, i, j, "GetMatEl");
nth(nth(cdr M, i), j) % element j of row i

; % GetMatEl

symbolic operator SetMatEl;
symbolic procedure SetMatEl(M, i, j, value);
% The assignment M := SetMatEl(M, i, j, value) effects
% M(i,j) := value where M is a LOCAL matrix.
% Based on setmatelem in matr.red
<<
LocMatCheck (M , i, j, "SetMatEl") ;
rplaca(pnth(nth (cdr M, i), j), value);
M

; % SetMatEl

symbolic operator LocMatCheck;
symbolic procedure LocMatCheck(M, i, j, Proc);
if not eqcar(M,'mat) then rederr list
(Proc, ": first argument not a matrix")
else if not fixp i then rederr list
(Proc, ":", i, "invalid as first matrix index")
else if not fixp j then rederr list
(Proc, ":", j, "invalid as second matrix index");

procedure Norm1(M);
% Returns 1-norm of matrix M =
% max. column sum of absolute values.
begin scalar norm, colsum; integer n ;
norm := length M; % list {nrows, ncols}
if (n := first norm) neq second norm then
rederr "Norm1: matrix not square" ;
norm := 0;
for j := 1 : n do <<
colsum :=
for i := 1 : n sum abs GetMatEl(M, i, j);
if colsum > norm then norm := colsum

;
return norm
end;

m1:=mat((1,2),(3,4));
norm1(m1); % ok
m2:=mat((sin(1),2),(3,4));
norm1(m2); % gives error

Does anyone know why the second norm calculation gives an error and how to fix it?

Thanks.

[Thomas Sturm]

Please, post your code within code tags, or attach a file.

Best,

Thomas

[miktexuser]

Here is the code thanks.

off int;

load_package matrix;

symbolic operator GetMatEl;
symbolic procedure GetMatEl(M, i, j);
% Returns M(i, j) where M is a LOCAL matrix.
% Based on getmatelem in matr.red
% Internally , M =(mat (m11 m12 ... )(m21 m22 ... ) ... )
<<
LocMatCheck(M, i, j, "GetMatEl");
nth(nth(cdr M, i), j) % element j of row i
>> ; % GetMatEl

symbolic operator SetMatEl;
symbolic procedure SetMatEl(M, i, j, value);
% The assignment M := SetMatEl(M, i, j, value) effects
% M(i,j) := value where M is a LOCAL matrix.
% Based on setmatelem in matr.red
<<
LocMatCheck (M , i, j, "SetMatEl") ;
rplaca(pnth(nth (cdr M, i), j), value);
M
>>; % SetMatEl

symbolic operator LocMatCheck;
symbolic procedure LocMatCheck(M, i, j, Proc);
if not eqcar(M,'mat) then rederr list
  (Proc, ": first argument not a matrix")
else if not fixp i then rederr list
  (Proc, ":", i, "invalid as first matrix index")
else if not fixp j then rederr list
  (Proc, ":", j, "invalid as second matrix index");

procedure Norm1(M);
% Returns 1-norm of matrix M =
% max. column sum of absolute values.
begin scalar norm, colsum; integer n ;
  norm := length M; % list {nrows, ncols}
if (n := first norm) neq second norm then
    rederr "Norm1: matrix not square" ;
norm := 0;
for j := 1 : n do <<
    colsum :=
       for i := 1 : n sum abs GetMatEl(M, i, j);
    if colsum > norm then norm := colsum
>>;
return norm
end;

m1:=mat((1,2),(3,4));
norm1(m1); % ok

m2:=mat((sin(1),2),(3,4));
norm1(m2); % gives error

[Eberhard Schruefer]

You haven't told us what your error looks like. But I assume it is something
like "* abs(sin(1)) + 3 invalid as number". You need to have
on numval,rounded; for the comparison in your procedure norm1 to work.

Or are refering to different error?

BTW, I changed REDUCE about a year ago to allow for local matrices and
matrix valued procedures. So in the newer version you could write

m := mat((sin(1),2),(3,4));

procedure Norm1(M);
% Returns 1-norm of matrix M =
% max. column sum of absolute values.
begin scalar norm, colsum; integer n ;
norm := length M; % list {nrows, ncols}
if (n := first norm) neq second norm then
rederr "Norm1: matrix not square" ;
norm := 0;
for j := 1 : n do <<
colsum :=
for i := 1 : n sum abs m(i,j);
if colsum > norm then norm := colsum>>;
return norm
end;

20: norm1 m;

* abs(sin(1)) + 3 invalid as number

21: on numval,rounded;

22: norm1 m;

6

Hope this helped.

Eberhard

[miktexuser]

Eberhard

You have fixed my problem.

I didn't know about your changes for local matrices. Does that mean that both getmatel and setmatel from MacCallum are not needed anymore?

If there a document that describes any changes such as this one that have been made to REDUCE recently?

Thanks.

[Eberhard Schruefer]

On 11/04/2015 12:38 AM, miktexuser wrote:

Eberhard

You have fixed my problem.

I didn't know about your changes for local matrices. Does that mean that both getmatel and setmatel from MacCallum are not needed anymore?

The functions getmatel/setmatel you refer to shouldn't be needed.

If there a document that describes any changes such as this one that have been made to REDUCE recently?

Unfortunately not. Below is my announcement to the developers list.

Eberhard

Subject: matrix valued procedures
Date: Sun, 07 Sep 2014 00:23:41 +0200
From: Eberhard Schruefer eschruefer@ca-musings.de
To: reduce-algebra-developers@lists.sourceforge.net

Prompted by a recent post to reduce-algebra:discussion I committed some rough conservative changes
to allow for matrix valued procedures (which is similar to what I did for lists). Unfortunately, due to its age
the matrix packages does not use the 'evfn' mechanism and things might be a little bit obscure.
Anyway, I put the code out there just to play with.

Matrix valued procedure are defined by the keyword 'matrixproc' (well 'matrix procedure' would have
been nicer, but I didn't want to complicate things).

Eberhard

Here are some examples

load matrix;

%Examples from the post to reduce-algebra discussion.

matrixproc SkewSym1 (w);
  mat((0,-w(3,1),w(2,1)), (w(3,1),0,-w(1,1)), (-w(2,1), w(1,1), 0));

skewsym1

X := SkewSym1(mat((qx),(qy),(qz)));

     [  0     - qz   qy  ]
     [                   ]
x := [ qz      0     - qx]
     [                   ]
     [ - qy   qx      0  ]

X * mat((rx),(ry),(rz));

[ qy*rz - qz*ry  ]
[                ]
[ - qx*rz + qz*rx]
[                ]
[ qx*ry - qy*rx  ]

SkewSym1(mat((qx),(qy),(qz))) * mat((rx),(ry),(rz));

[ qy*rz - qz*ry  ]
[                ]
[ - qx*rz + qz*rx]
[                ]
[ qx*ry - qy*rx  ]

matrixproc quatInv4(a,b,c,d);
  mat((a),(-b),(-c),(-d));

quatinv4

matrixproc quatInv1(p);
  quatInv4(p(1,1),p(2,1),p(3,1),p(4,1));

quatinv1

q:=mat((qs), (qx), (qy), (qz));

     [qs]
     [  ]
     [qx]
q := [  ]
     [qy]
     [  ]
     [qz]

quatInv4(q(1,1),q(2,1),q(3,1),q(4,1));

[ qs  ]
[     ]
[ - qx]
[     ]
[ - qy]
[     ]
[ - qz]

quatInv1(q);

[ qs  ]
[     ]
[ - qx]
[     ]
[ - qy]
[     ]
[ - qz]

% Some mor examples.

matrixproc foo u;
  begin
    u(1,1) := 1;
    u(2,1) := a*u(3,1) + u(2,1);
    return u
  end;

foo

foo q + q;

[  qs + 1   ]
[           ]
[a*qy + 2*qx]
[           ]
[   2*qy    ]
[           ]
[   2*qz    ]

m := mat((a,b),(c,d));

     [a  b]
m := [    ]
     [c  d]

id := mat((1,0),(0,1));

      [1  0]
id := [    ]
      [0  1]

matrixproc exmpl2 u;
  u**2 - trace u*u;

exmpl2

exmpl2 m/det m;

[-1  0 ]
[      ]
[0   -1]

[miktexuser]

Thank you Eberhard. Matrixproc is very useful.