@Markiyan Hirnyk 

I did not read them because as I have added (probably you did not see because it was during your reply) I am interested in general Groebner bases which are inherently huge.

Edit. OK, you made me curious and I have looked over [1]. As I have anticipated, they approximate a Groebner basis using Newton's method. That is not what most people want!

@Markiyan Hirnyk 

OK, but if for 70 equations we need supercomputers then ...

Edit. I was (and am) interested in general Groebner bases, not in some versions designed for numerics. Here the situation is unfortunately clear: the Grobner basis is unique and very often is HUGE.

@Markiyan Hirnyk 

The Groebner basis technique is excellent in theory and works for small size problems. I like it very much.
Unfortunately for even moderately large problems it cannot be used. The compexity could be double exponential O(a^(b^n))   !!

@Carl Love 

Yes, irem is a little faster than mods, but if p and nops(S) are huge,
taking (at least from time to time) the irem of the partial product
will speed up the computations.

@acer 
Probably a more instructive example would be:

progbar:= proc(p::realcons)
local x;
if p>100 or p<0 then DocumentTools:-Tabulate( [[]] ): return NULL fi:
DocumentTools:-Tabulate([ sprintf("%a %c", p, "%"),
                          plot(100,x=0..p, view=[0..100,0..100],
                                  filled,axes=box, tickmarks=[[],[]],labels=["",""], size=[450,50]) ],
                          exterior=none, interior=none, widthmode=pixels, width=450 );
end:

s:=0;
"starting ...";
for i to 10000 do
  s:=s + evalf(1/i);
  if i mod 100 = 0 then  progbar(i/100) fi;
od:
progbar(101):
's'=s;

@dpaddy 

If g is continuous and C^1 piecewise continuous then you can use IntegrationTools[Split].
If g is given as piecewise, Maple should be able do compute directly int(f*D(g), a..b);

If g is C^1 piecewise continuous but not (globally) continuous you must also consider the jumps.

@dpaddy 

If in the integral of f*dg  one of the functions  f,g  is C^1 (at least piecewice) then the integral can be reduced to a Riemann one
and it is not difficult to write a Maple program to do this.
(But be aware of the fact that there exists (pathological) functions for which Maple fails in computing the Riemann integral.)

Otherwise, you have to study theoretically the problem, and if the integral exists try to reduce it to the previous case (if possible).
Note that the Stieltjes integral is not very simple in the general case. Even if both f,g are continuous, it may not exist.

@AmusingYeti 

Any 3d plot!

plot3d(x^2-y^2,x=-1..1, y=-1..1);
#Here is the print screen; I have posted such images here long time ago.
(I have to  use Maxima for vector 3d plots!)

@AmusingYeti 

It must be added that unfortunately the .eps export for 3d graphics is terrible (unusable), at least in Windows.

@Carl Love 

g(x)=floor(x):

inf(x*dg(x), x=0..5) =

x*g(x)|(x=5) - x*g(x)|(x=0)  - int(g(x), x=0..5) =

25 - 10 =

15

Now, using Maple:

Int(x*D(g)(x), x=0..5):
IntegrationTools[Parts](%,x):
eval(%, g=floor):
value(%);
                               15

# Wrong result if trying directly!
int(x*D(floor)(x), x=0..5);
                               0

@Carl Love 

The Stieltjes integral  Int_[0,5]  x d(floor(x))    cannot be computed by Maple. It must be transformed by hand.

@Carl Love 

Thank you for the interesting info.

Best regards,
V.

@Joe Riel 

Thank you, I found it.
Unfortunately it seems that for an interested user it will be harder and harder to know what is Maple actually doing. It used to be so nice to see the Maple code when something went wrong or, on the contrary, suprisingly well!

@Tyna 

V is supposed to be a vector.

V^+  converts this vector from row to column (or vice-versa).

A and B were filled with some values (by RandomMatrix) in order to obtain numeric results
(you did not provide the numeric matrices).

If V is a matrix (instead of a vector), the objective function is not a scalar! You cannot divide two matrices or minimize a matrix!
[Actually this would be possible in vector optimization, but then a lot of details must be considered.]

So, you should post the whole problem (or a smaller version) if you want a numeric solution.
But it has to make sense :-)

@Joe Riel 

Yes, I forgot about _rest. I have used myself initially args[3..-1] but it was slower.

(But actually the timing seems to be rather inconsistent here.)