Speed up bpy object creation with "for loop" so i'm fairly new to blender scripting and could need some advice with this:

To put it simply, I want to create a cube made up of smaller cubes (100x100x100) with this script:

N=100;
for i in range(0,N,3):   
   for j in range(0,N,3):   
      for k in range(0,N,3):
          bpy.ops.mesh.primitive_cube_add(enter_editmode=True, location=(i,j,k))   


which works well with smaller N. But takes horribly long with larger N. In the code above there are over 9000 objects to be created so I can assume there is a lot of memory to be allocated which slows things down.

However, I there a way to create such objects much faster? Is there a fundamental problem with my approach that I've not seen yet?

I also experience this "slow" behavior in other scripting modellers like OpensScad and hoped maybe blender does this a little better.

Tips and advice very much appreciated! Thanks!
Using array modifier

Using array modifiers will speed up the processing because it will lower the overall io between Blender and Python and also lower bpy.ops usage:

Create a base cube object.

Create and apply three array modifiers on it.

The code is commented below, but if not clear enough, please ask in the comments.

import bpy
import bmesh

def apply_modifier( obj, count, vector ):
    # Add an array modifier with its wanted parameters
    m = obj.modifiers.new( name = "Array", type = 'ARRAY' )
    m.use_relative_offset = False
    m.use_constant_offset = True
    m.constant_offset_displace = vector
    m.count = count
    # Apply it
    bpy.ops.object.modifier_apply( modifier = m.name )

def create_base_cube( context, name, size ):
    # Create an object
    mesh = bpy.data.meshes.new( name )
    obj = bpy.data.objects.new( name, mesh )

    # Add a cube
    bm = bmesh.new()
    bm.from_mesh( mesh )
    bmesh.ops.create_cube( bm, size = size, calc_uvs = False )
    bm.to_mesh( mesh )
    bm.free()

    # Link the object to the scene and set it active
    context.scene.collection.objects.link( obj )
    context.view_layer.objects.active = obj

    return obj

name = "CubeMap"
size = 2

obj = create_base_cube( bpy.context, name, size )

N=100;
interval = 3

apply_modifier( obj, N, (interval, 0, 0) )
apply_modifier( obj, N, (0, interval, 0) )
apply_modifier( obj, N, (0, 0, interval) )


Using dupliverts (proposed by Jackdaw in the comments)

An alternative is to use dupliverts, which seems to be faster and allows to define any patterns (not only cubic).

The principle is to create an object with vertices at the positions you want, then make it parent to a cube and use dupliverts as instancing method.

import bpy
import bmesh
from itertools import product

def create_base_cube( context, name, size ):
    # Create an object
    mesh = bpy.data.meshes.new( name )
    obj = bpy.data.objects.new( name, mesh )

    # Add a cube
    bm = bmesh.new()
    bm.from_mesh( mesh )
    bmesh.ops.create_cube( bm, size = size, calc_uvs = False )
    bm.to_mesh( mesh )
    bm.free()

    # Link the object to the scene and set it active
    context.scene.collection.objects.link( obj )
    context.view_layer.objects.active = obj

    return obj

def make_result_object( context, name, cube, vertices ):
    mesh = bpy.data.meshes.new( name )
    obj = bpy.data.objects.new( name, mesh )

    # Associate the vertices
    obj.data.from_pydata( vertices, [], [] )

    # Link the object to the scene and set it active
    context.scene.collection.objects.link( obj )
    context.view_layer.objects.active = obj

    # Make the object parent of the cube
    cube.parent = obj
    # Make the object dupliverts
    obj.instance_type = 'VERTS'

    return obj

name = "CubeMap"
size = 2

cube = create_base_cube( bpy.context, name, size )

N=3*100
interval = 3

vertices = [(x, y, z) for x, y, z in product( range(0,N,interval), range(0,N,interval), range(0,N,interval) )]

result = make_result_object( bpy.context, "result", cube, vertices )