Assign a vector using another vector

Asked by James Maddison

I would like to assign a vector to be equal to another vector. Something like:

F1.set_local(F2.array() + F3.array())
F1.apply(...)

is quite expensive compared to a simple:

F1 = F2 + F3

How do I assign a (previously created) vector F1 to be equal to the sum of F2 + F3 without using the array method? I've profiled the latter and found this to be many times slower than the direct addition of two vectors, presumably due to extra copying. Would it be possible for a method of the form:

F1.assign(F2 + F3)

to be added to vectors?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
Last query:
Last reply:
Revision history for this message
Best Johan Hake (johan-hake) said :
#1

On 09/19/2012 04:51 PM, James Maddison wrote:
> New question #209010 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/209010
>
> I would like to assign a vector to be equal to another vector.
> Something like:
>
> F1.set_local(F2.array() + F3.array()) F1.apply(...)
>
> is quite expensive compared to a simple:
>
> F1 = F2 + F3
>
> How do I assign a (previously created) vector F1 to be equal to the
> sum of F2 + F3 without using the array method? I've profiled the
> latter and found this to be many times slower than the direct
> addition of two vectors, presumably due to extra copying. Would it be
> possible for a method of the form:
>
> F1.assign(F2 + F3)
>
> to be added to vectors?

You should be able to do this already:

  F1[:] = F2 + F3

as F1[:] maps to the assignment operator.

Or even better:

  F1[:] = F2
  F1 += F3

Johan

Revision history for this message
James Maddison (jamesmadd) said :
#2

That does indeed seem to work, although the assignment is surprisingly slow. e.g.:

#!/usr/bin/env python

from dolfin import *
import time

mesh = Rectangle(0.0, 0.0, 1.0, 1.0, 1000, 1000)
space = FunctionSpace(mesh, "CG", 1)
F1 = Function(space); F1.assign(Constant(1.0))
F2 = Function(space); F2.assign(Constant(1.0))

cl = time.time()
for i in range(1000):
  F3 = F1.vector() + F2.vector()
print "%.6f" % (time.time() - cl)

cl = time.time()
for i in range(1000):
  arr = F1.vector().array() + F2.vector().array()
  F3.set_local(arr)
  F3.apply("")
print "%.6f" % (time.time() - cl)

cl = time.time()
for i in range(1000):
  F3[:] = F1.vector() + F2.vector()
print "%.6f" % (time.time() - cl)

cl = time.time()
for i in range(1000):
  F3[:] = F1.vector()
  F3[:] += F2.vector()
print "%.6f" % (time.time() - cl)

yields:

3.332180
19.223399
4.431610
5.573048

The assignment (third case) adds an extra 30% to the runtime (over the first case).

Revision history for this message
James Maddison (jamesmadd) said :
#3

Thanks Johan Hake, that solved my question.

Revision history for this message
Johan Hake (johan-hake) said :
#4

On 09/19/2012 06:01 PM, James Maddison wrote:
> Question #209010 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/209010
>
> James Maddison posted a new comment:
> That does indeed seem to work, although the assignment is surprisingly
> slow. e.g.:
>
> #!/usr/bin/env python
>
> from dolfin import *
> import time
>
> mesh = Rectangle(0.0, 0.0, 1.0, 1.0, 1000, 1000)
> space = FunctionSpace(mesh, "CG", 1)
> F1 = Function(space); F1.assign(Constant(1.0))
> F2 = Function(space); F2.assign(Constant(1.0))
>
> cl = time.time()
> for i in range(1000):
> F3 = F1.vector() + F2.vector()
> print "%.6f" % (time.time() - cl)
>
> cl = time.time()
> for i in range(1000):
> arr = F1.vector().array() + F2.vector().array()
> F3.set_local(arr)
> F3.apply("")
> print "%.6f" % (time.time() - cl)
>
> cl = time.time()
> for i in range(1000):
> F3[:] = F1.vector() + F2.vector()
> print "%.6f" % (time.time() - cl)
>
> cl = time.time()
> for i in range(1000):
> F3[:] = F1.vector()
> F3[:] += F2.vector()
> print "%.6f" % (time.time() - cl)
>
> yields:
>
> 3.332180
> 19.223399
> 4.431610
> 5.573048
>
> The assignment (third case) adds an extra 30% to the runtime (over the
> first case).

Try:

  for i in range(1000):
    F3[:] = F1.vector()
    F3 += F2.vector()
  print "%.6f" % (time.time() - cl)

By using in place addition to a slice (F3[:]+=F2.vector()) you invoke
some peculiar Python magic. What you wrote in the last for loop equals to:

for i in range(N):
  F3[:] = F1.vector()
  F4 = F3[:] # Which triggers a copy!
  F4 += F2.vector()
  F3[:] = F4

Which comes with a performance hit over my suggestions.

Johan

Revision history for this message
James Maddison (jamesmadd) said :
#5

Yes, you're right, your suggestion yields a time of 2.725569s.